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ABSTRACT 

We present the first three-dimensional magnetohydrodynamic (MHD) simu- 
lations of a circumbinary disk surrounding an equal mass binary. The binary 
maintains a fixed circular orbit of separation a. As in previous hydrodynamical 
simulations, strong torques by the binary can maintain a gap of radius ~ 2a. 
Streams curve inward from r ~ 2a toward the binary; some of their mass passes 
through the inner boundary while the remainder swings back out to the disk. 
However, we also find that near its inner edge the disk develops both a strong 
m — 1 asymmetry and growing orbital eccentricity. Because the MHD stresses 
introduce more matter into the gap, the total torque per unit disk mass is ~ 14 
times larger than found previously. The inner boundary accretion rate per unit 
disk mass is ~ 40 times greater than found from previous hydrodynamical calcu- 
lations. The implied binary shrinkage rate is determined by a balance between 
the rate at which the binary gains angular momentum by accretion and loses it 
by gravitational torque. The large accretion rate brings these two rates nearly 
into balance, but in net, we find that a/ a < 0, and its magnitude is about 2.7 
times larger than predicted by the earlier hydrodynamic simulations. If the bi- 
nary comprises two massive black holes, the accretion rate may be great enough 
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for one or both to be AGN, with consequences for the physical state of the gas 
both in the disk body and in its inner gap. 

Subject headings: accretion, accretion disks - - binaries: general — MHD - 
methods: numerical 



1. INTRODUCTION 

Various types of astronomical binary systems can be embedded in gaseous disks, from 
young binary stars to stars with growing planets to binary black holes. Such disks have been 
observed directly in nearby star-forming regions. One of the best resolved (around a young 
binary) is in GG Tau (e.g., Dutrey et al. 1994; Krist et al. 2005). To date, however, only a 
few possible disk-planet systems have been directly imaged (e.g., Greaves et al. 2008; Kalas 
et al. 2008; Hashimoto et al. 2011). Although there is no direct evidence for the existence 
of circumbinary disks involving binary black hole systems, it is generally believed that such 
configurations should exist near the centers of galaxies after a galaxy merger (e.g., Begelman 
et al. 1980; Ivanov et al. 1999; Merritt & Milosavjevic 2005; Escala et al. 2004, 2005; 
Mayer et al. 2007; Dotti et al. 2009). 

Tidal forces exerted by the binary can sometimes clear a gap in the disk. When the 
binary mass ratio q = M2/M1 <C 1, where Mi denotes the mass of the star and M2 the 
mass of the planet, the gap that is formed is an annular ring around the primary through 
which the secondary travels. Whether such a gap opens depends on whether the secondary's 
mass is sufficiently large to overcome the gap closing effects of internal stresses. Independent 
of whether such a gap exists, the secondary can exchange angular momentum with the gas 
via gravitational torques. Inward orbital migration of the secondary may occur on the disk 
inflow timescale if the disk mass is large compared to the secondary mass (Lin & Papaloizou 
1986, 1993). Otherwise the migration will be slower (Ivanov et al. 1999; Armitage & 
Natarajan 2002). There has been extensive theoretical study of this situation, using both 
analytic and numerical methods (e.g., Goldreich & Tremaine 1980; Lin & Papaloizou 1986, 
1993; Bryden et al. 1999; Ivanov et al. 1999; Bate et al. 2003; Nelson & Papaloizou 2003; 
Winters et al. 2003). 

When q is closer to unity, the gap can include the entire binary itself. In this case, the 
resulting configuration can contain as many as three disks: one around each member of the 
binary and one that orbits outside the binary, called the circumbinary disk. Observational 
evidence of such large gap clearing and circumbinary disks has been found in several young 
stellar binaries (Mathieu 1994). Numerical simulations have also been applied to study 
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q < 1 binary systems with a circumbinary disk (e.g., Artymowicz & Lubow 1994; Bate 
k Bonnel 1997; Giinther & Kley 2002; Escala et al. 2005; MacFadyen & Milosavljevic 
2008; Hayasaki et al. 2007; Cuadra et al. 2009; Hanawa et al. 2010). These studies have 
found that the radius of the disk inner edge depends on several factors, including the binary 
separation, mass ratio, eccentricity, and the strength of the disk turbulence (Artymowicz & 
Lubow 1994). 

In a one-dimensional model (one that considers only radial dependences), a circumbi- 
nary disk would in principle behave as a "decretion" disk' (Pringle 1991): the binary loses 
angular momentum to the surrounding disk so that the disk is repelled. However, the one- 
dimensional model neglects any non-axisymmetric properties of the disk. For instance, there 
could be some directions where the binary torque becomes weak enough that matter could 
leak through the disk inner edge. The gas then penetrates the gap and may accrete onto 
the binary. Indeed, this penetration effect has been witnessed in numerous two- or three- 
dimensional simulations of disks with various binary mass ratios using either smoothed par- 
ticle hydrodynamics (SPH; Artymowicz & Lubow 1996; Escala et al. 2005; Cuadra et al. 
2009; Dotti et al. 2009) or grid-based methods (Bryden et al. 1999; Giinther & Kley 2002; 
MacFadyen & Milosavljevic 2008; de Val-Borro et al. 2011). Such work has found that 
in the low-density gap there are gas flows from the disk inner edge to the binary in form of 
narrow, high-velocity spiral streams. The flow rate through the gap depends on the binary 
and the disk properties. Compared to the accretion rate near the disk center that would 
be expected in the absence of binary torque, the accretion rate appears to be reduced (e.g., 
Lubow et al. 1999; Lubow & D'Angelo 2006; MacFadyen & Milosavljevic 2008, hereafter 
MM08). However, Rozyczka & Laughlin (1997) found no reduction at all. 

Disk eccentricity is another potential non-axisymmetric property. A disk can become 
eccentric as a result of its interaction with the binary. Circumbinary disks can, of course, 
gain eccentricity through direct driving by an eccentric binary (Artymowicz et al. 1991; 
Papaloizou et al. 2001; Rodig et al. 2011). However, there is also evidence that disk 
eccentricity can arise even when the orbit of the binary is circular (e.g., Papaloizou et al. 
2001, MM08). Simulations of protoplanetary disks have shown that these disks are subject to 
a resonant mode coupling instability (Kley & Dirksen 2006; D'Angelo et al. 2006). Through 
this instability, the disk's eccentricity can grow although the planet is on a fixed circular 
orbit. This instability follows the same tidal resonance mechanism found for eccentrically 
unstable circumstellar disks in superhump binaries (Lubow 1991; Kley et al. 2008). For 
circumstellar disks, Lubow (1994) found that orbiting secondaries can drive eccentricity by 
stream impacts if they are strongly modulated in time. Recently, MM08 found that the disk 
around an equal mass binary on a circular orbit also became eccentric after a large number of 
binary orbits. They suggested an eccentricity generation mechanism that involved the action 
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of the binary torque on the gas within the low density gap. Conversely, disk eccentricity might 
also excite eccentricity in the binary (e.g., Papaloizou et al. 2001). For coalescing massive 
black holes, the residue eccentricity in the emitted gravitational waves might be detected by 
proposed instruments like the Laser Interferometer Space Antenna (LISA), perhaps signaling 
gas-driven evolution (Armitage, & Natarajan 2005; Key & Cornish 2011). 

Almost all previous studies of circumbinary disks have adopted the a-prescription to 
describe internal stresses and treat either as evolving through diffusion (in one dimension) 
or as a result of "viscous" stresses (in two or three dimensions). In these efforts, the internal 
stress to pressure ratio a was taken to be constant everywhere and at all times. Although that 
might be a reasonable approximation for vertically-integrated and time-averaged conditions 
in the main body of a disk, it becomes unrealistic for highly time-dependent turbulent 
accretion flows and low density regions outside the disk body(Hawley & Krolik 2001). 
Since the exchange of angular momentum between the binary and the disk is crucial for 
both the circumbinary disk and the binary, we need more a realistic description of the 
underlying internal stresses. It is now generally recognized that whenever the material of 
the disk is sufficiently ionized so as to be well-coupled to any embedded magnetic field, 
the principal mechanism of angular momentum transport is MHD turbulence induced by 
the magnetorotational instability (MRI). It is therefore necessary to study circumbinary 
disks using MHD simulations in which internal stress arises self-consistently from turbulence 
generated by the MRI. To date, this has been done only for extreme mass ratio star-planet 
systems, using either unstratified (Nelson Sz Papaloizou 2003; Winters et al. 2003; Baruteau 
et al. 2011) or stratified (Uribe et al. 2011) MHD simulations. 

It is the goal of this project to construct the first three-dimensional (3D) MHD simula- 
tion of a circumbinary disk around an equal mass binary. To simplify our model, we assume 
the disk and binary to be coplanar. The binary orbits on a fixed circular orbit. On the basis 
of this simulation, we will try to answer the following questions: (1) What is the inner disk 
structure? Is the disk truncated or not? (2) How is angular momentum transported within 
the disk? Can the internal stress balance the binary torque and therefore allow the disk to 
achieve a quasi-steady state? (3) Is there any eccentricity growth of the disk? If so, what is 
the cause? (4) How does the accretion rate onto a binary compare with the rate onto single 
point mass? 

We organize this paper as follows: In § 2, we describe the physical model and numerical 
procedures of our circumbinary disk simulations. In § 3, we present our simulation results. 
We then discuss the binary torque and binary contraction in § 4.1 and 4.2. In § 4.3, we 
discuss possible mechanisms for disk eccentricity growth. We explain the formation of an 
asymmetric density concentration near the disk inner edge in § 4.4. Finally, in § 5 we 
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summarize our conclusions. 



2. NUMERICAL SIMULATION 

In this section, we discuss in detail the numerical procedure of this work. The code used 
is a modern version of the 3D, time-explicit Eulerian finite-differencing ZEUS code for MHD 
(Stone & Norman 1992a,b; Hawley & Stone 1995). We modified the code to cope with the 
time- dependent binary potential. 

2.1. Physical Model 

We construct our physical disk model in the inertial frame in which the center of mass 
of the binary is at rest at the origin. Treating the simplest case first, we assume an equal 
mass circular restricted binary system, i.e. the binary orbits circularly in the disk plane, and 
we neglect binary evolution as the disk inflow time scale is usually much smaller than the 
shrinkage timescale of the binary, and certainly more than the duration of the simulation 
(see estimate in § 5). We set the gravitational constant G, total binary mass M and the 
binary separation a to be unity, and therefore the binary frequency f2bin = \j GM/a 3 is 
unity as well. We also assume the circumbinary disk to be cold and thin. As we are 
mainly concerned with the orbital dynamics of the flow, we choose a simple global isothermal 
sound speed c s = 0.05. The disk flares at larger radii because the ratio of height to radius 
H/R = Cs/RQk = 0.05(i?/a) 1 / 2 , where H denotes the scale height of the disk, R is the 
disk radius in cylindrical coordinates, and f2x = y GMJ R 3 is the Keplerian frequency. 
As the fluid is well coupled to the magnetic field even in a cold disk where the ionization 
fraction is far below unity(see, e.g., Stone et al. 2000), we assume ideal MHD. For this 
reason, we include no explicit diffusivity except the von Neumann-Richtmyer bulk viscosity 
in compressive regions that ensures the right jump conditions for shock waves. The effect 
on damping the angular momentum is negligible compared to other transport mechanisms. 
Because we are most concerned with the inner part of the disk, not accretion onto the 
binary, we excise a central region. This cut-out must be well inside the inner edge of the 
disk so that it does not affect fully resolving the inner disk and any gas leakage from the 
inner edge. We choose to cut out the area within 0.8a. This region is beyond the main 
extent of the interior disks that surround each binary member because these disks are each 
tidally truncated at ~ 0.3a from their central objects (Paczyhski 1977). Once the disk 
reaches the quasi-steady stage, we find that the inner edge of the disk is located at r ~ 2a, 
which is about a factor of 2.5 outside the cut-out. We approximate the time-dependent 
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binary potential in the disk region using Newtonian dynamics because the disk region we are 
interested in is far from the gravitational radii of the individual binary components. The disk 
mass is assumed to be much smaller than the binary mass, so that the Toomre parameter 
satisfies Q ~ (H/R){M/M$) 3> 1, where Md denotes the disk mass. We therefore neglect 
the contribution of the disk self-gravity to the potential. 



The properties of circumbinary disks require us to resolve three different length scales 
when setting up the grid. The first is the half disk thickness H. The computational domain 
has to contain at least several (< 4) scale heights on each side of the midplane, and for each 
H, several tens of cells are needed. The second is the maximum growth-rate wavelength 
of the MRI, Amri = 8n/\/l5v a/£1(R), where va is the Alfven speed, and Cl(R) is the disk 
rotational frequency. We require Amri to be resolved by at least six grid elements. The last 
one is the spiral density wavelength, Ad ~ 27rc s /flbm- There should be many cells across this 
wavelength (MM08); we require at least six. 

In order to satisfy the above requirements, we follow the scheme proposed in Noble et 
al. (2010) to construct our grid in spherical coordinates (r,9,<fi). We adopt a logarithmic 
grid in the radial direction, which provides a constant Ar/r. The vertical grid is derived by 
mapping a simple linear function y(x) = x for x G [0, 1] through a polynomial transformation 
(see equation (6) in Noble et al. (2010)): 



where £, 6 C and n are parameters that define the shape of the polynomial. Note that y = 0.5 
is exactly mapped to the midplane. The merit of this #-grid is that for n > 1 it ensures 
dense and nearly uniform cell elements close to the disk midplane. The azimuthal grid is 
evenly spaced and covers all 2tt. Theoretically speaking, the aspect ratios of the cell shape 
should be as isotropic as possible, but the azimuthal cell widths can be a factor of a few 
longer than the other two, as the shearing tends to draw out features in this direction. 

The grid resolution used in the present simulation is 400 x 160 x 540 in (r,9,(f)), with 
a computational domain covering [r in ,r out ] radially, [8 c ,tt — 8 C ] meridionally and [0, 27r] az- 
imuthally (see also in Table 1), where r in = 0.8a, r out = 16a, and 8 C = 0.2 is the width of 
the cut-out around the polar axis. The other parameters used in equation (1) are £ = 0.9 
and n — 9. Using this grid, we are able to resolve one disk scale height with 20 cells at the 
inner boundary of the computational domain. The resolution grows to ~ 40 cells per scale 



2.2. Grid Scheme and Boundary Conditions 
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height at radius ~ 3a. Within two scale heights of the midplane, Amri is resolved by more 
than six cells if the plasma /3 is no greater than ~ 100. The grid also resolves Ad with at 
least six cells in the r-0 plane for r < 7a. 

We choose outflow boundary conditions in the radial and meridional directions for the 
gas, which permit only flows going outward; any inward velocities are set to zero. As we 
cover all 2tt of 0, the boundary conditions for 0-grid are simply periodic. For the magnetic 
field in the radial and meridional directions, we set the transverse components of the field 
to be zero in the ghost zones. The components normal to the boundaries are calculated by 
imposing the divergence- free constraint. 



We begin with a prograde disk orbiting in the binary plane. Between r min = 3a and 
r maj( = 6a, its density is constant in the midplane. The initial disk is axisymmetric, and 
the polar angle density distribution is p = po ex P [ — {& — 71 /2) 2 / (\/2H /r) 2 ], in which p — 1 
is the unit of disk density. This form provides initial hydrostatic balance vertically for a 
point mass potential and zero radial pressure gradient along the midplane. For a first order 
approximation, the difference between a binary potential and a point-mass in the midplane 
can be described by the temporally and azimuthally averaged quadrupole moment of the 
binary potential. We therefore modify the angular frequency of the initial disk to account 
for the quadrupole contribution: 



where q = 1 is the mass ratio of the binary. Here we replace R with r as r = R sin 6 R for 
regions near the midplane. We also add l/(rH)dP/dr to the right hand side of equation (2) 
to compensate for the small radial gradient of the vertically integrated pressure P = £c^. 

The initial magnetic field is a single poloidal loop within the main body of the disk. 
It is subthermal with plasma /3 = 100 on average. The field is constructed from the vector 
potential A = (0,0, A^) and we define A^ by A^ = A 0y /psm(2iir/kH)(r/r min — 1)(1 — 
r/r mM ) — y/p cu t if At/, > and zero otherwise, where k = 2Q hin a/c s , p cut = 10 _3 po, and A is 
a constant determined by the constraint on the averaged (3. 



2.3. Initial Setup 




(2) 
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2.4. Diagnostics 

Three-dimensional MHD simulations usually produce a large amount of data, which 
poses challenges for data storage, transport and post-simulation analysis. In order to fa- 
cilitate the study of the spatial and temporal properties of the disk, we write out spatially 
averaged history data at a frequency of one dump per time unit and write out 3D snapshots 
every five units of time. 

We use two different types of spatially averaged history data. The first is defined by 
either an integration or an average over radial shells. Shell-averaging for variable X is defined 
by 

. . „ f Xr 2 sin OdOdd) , , 



where J r 2 sin 8d8d<p = A^(r) is the shell surface area. With this definition, the density- 
weighted shell average is (X) p = (pX)/(p). For example, the net disk accretion rate is 



M(r,t) = J pv r r 2 sin 6d0d(j> = A x (pv r ), (4) 

and the average specific angular momentum / = (^r) p = (pv^r) / (p). The surface density is 
vertically integrated and azimuthally averaged quantity: 

E(r,t) = — I pr sin 6d6d(j). (5) 
2it J 

The second type of average is two-dimensional, either an azimuthal average of poloidal 
slices or a vertical average referred to the equatorial plane. We define the vertical average 
by 

r . . , . f Yr sin 9d9 , , 

(YUr,<P,t) = J . , (6) 

J r smvdo 

and the density weighted vertical average by 

(YM, t ) = lf^ =^. (7) 
J prsmvaV {p) z 

Similarly we have the azimuthal average is (Y)^(r, 9,t) = ^- J Yd(j). 

We need to be very careful about the definition of the r-0 component of the internal 
stress in the present simulation. We follow Haw ley & Krolik (2001) and define the stress as 

B B 

w r Jr, 9, 0, t) = t r s + r r s = r — - + p8v r 8vs, (8) 

47T 
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where t r( p is the Maxwell stress and r r<i> is the Reynolds stress. The perturbed velocities 5v r 
and 5v<p are calculated from 



The vertically integrated and azimuthally averaged total stress can then be described by 



where T r( p and _R r </> are the average Maxwell stress and Reynolds stress respectively, and 
L z = A x /(2irr) is the vertical integral length, which in our case equals the height of the 
computational domain at a given radius r. In a similar way, we can obtain the vertically 
averaged stress (w r( p) z and azimuthally averaged stress (w^)^ as well. 



To study the mechanism for disk eccentricity growth, we also carried out a set of numer- 
ical experiments using two-dimensional viscous hydrodynamic simulations. Their purpose 
was both to distinguish hydrodynamic from MHD effects and to test the effect of where the 
inner boundary is placed. 

The hydrodynamic simulations in this paper used the a-disk prescription. We chose 
a — 0.1 throughout the disk because that is roughly the ratio of stress to pressure in the 
disk body of our MHD simulation (note MM08 used a = 0.01). The hydrodynamic disks are 
evolved in polar (r, 0) coordinates in the inertial frame. Following the grid scheme of the 
MHD simulation, we set the radial grid to be logarithmically spaced, while the azimuthal 
grid is spaced uniformly. 

Among these simulations, B2D.rin=0.8 serves as the control run. The parameters which 
describe the physical properties of the hydrodynamic disk and the binary are kept the same 
as in the MHD case. The resolution of the control run is 512 x 1024 for r x 0. Its grid 
covers the same physical extent in radial and azimuthal directions as the MHD one. The 
initial surface density of the two-dimensional disk is simply taken from a vertical integration 
of the initial condition of the MHD disk. Other hydrodynamic simulations are reruns of 



5v r (r, 9, 0, t) 
Sv^r, 9,<f),t) 



v r - (v r ) p 
V4, - {v^p. 



(9) 




(10) 



2.5. Hydrodynamic Simulations 
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B2D.rin=0.8 at t — 500 with various locations of the inner boundary. The initial disks of 
the reruns are obtained by truncating the restart data of B2D.rin=0.8 to the desired radius 
while keeping Ar/r and A0 fixed. The properties of both the reruns and the control run 
can be found in Table 1. 



3. RESULTS 

We present one 3D MHD simulation, called B3D, in this paper. 1 This simulation was 
terminated after the inner disk (r < 3a) reached a quasi-steady stage for several hundred 
time units. Longer simulations require better resolution to resolve the MRI in the growing 
density concentration region (see § 3.3.3). 

B3D ran for ~ 480 code units, which corresponds to ~ 77 binary orbits. The gas quickly 
fills in the initially empty region within several binary orbital periods, and after ~ 100 units 
of time, the disk becomes fully turbulent. The binary torque then is able to maintain a low 
density gap, and the inner part of the disk finally reaches a quasi-steady state after t ~ 200. 
We note that this single simulation consumed ~ 720K CPU hours on the Kraken Cray XT5 
system. 

In the following subsections, we first describe the overall evolution of the circumbinary 
disk and its quasi-steady state. We then try to discuss how the angular momentum is 
transported in the inner disk in § 3.2. We will mainly discuss the characteristic disk structures 
in § 3.3 and 3.4. The field structure will be considered in § 3.5 and the temporal properties 
of the accretion in the last subsection. 



3.1. A Secularly Evolving Quasi-steady State 

After the first 100 time units, the binary torque starts to clear out a low density gap 
between the inner boundary and r ~ 2a. In the top- left panel of Figure 1, we show the 
vertically integrated surface density of the circumbinary disk at t — 120 in the x — y plane. 
In that panel, we can clearly see the disk is truncated at around twice the binary separation. 
We also find two streams emanating from the disk edge toward the binary components. From 
the other three panels in Figure 1, we find that as the simulation continues the gap persists, 



1 We also performed a short duration rerun for t = 300-322 with higher dump rates: ten history 
dumps and one 3D dump per time unit. They are used when high time resolution is required, e.g. 
when we try to investigate the angular momentum budget and the stream dynamics. 
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the surface density gradually increases outside the gap, and finally an incomplete ring of 
dense gas forms due to the combined effects of mass accretion and gravitational torque. The 
stream also persists, but only one arm at a time. We will study the properties and effects of 
the transient stream in § 3.4 and 4. At t = 250 the disk appears to be elliptical and slightly 
off center, and at t — 350 and 450 the disk obviously becomes eccentric. We also observe a 
growing azimuthal density asymmetry. We will discuss the eccentricity growth and density 
asymmetry in § 3.3 and 4. Despite the slow growth of both disk eccentricity and asymmetry, 
the principal disk structures, the gap and stream, are sufficiently time-steady to suggest that 
in a qualitative sense the simulation has achieved a statistically stationary state. An analysis 
based on angle-averaged quantities will further demonstrate that. 

We show the enclosed mass (integrated disk mass interior to a given radius) inside r = a, 
1.5a, 2a, 3a and 4a as a function of time in Figure 2. We find that the mass at r < 3a (i.e., 
inside the initial disk's inner edge) undergoes dramatic growth during the first ~ 100 time 
units as the initial disk fills in regions with < r min . The radial pressure gradient at the edge 
of our initial disk leads to an inflow during the first several orbits, and once toroidal field 
develops near the disk inner edge, the Maxwell stress quickly removes the angular momentum 
of the low density flow and drives inflow. As this transient phase passes, the mass in the 
gap levels off and stays quasi-steady (falling very slowly) after t ~ 200. We therefore define 
the quasi-steady stage as the period between t = 200 and the end of the simulation. On the 
other hand, the mass inside r = 3a and 4a rises by a factor of two after t = 100, indicating 
that this region does not achieve inflow equilibrium. The mass-interior profiles also possess 
high frequency fluctuations (typical time scale ~ 3 time units) due to the orbital modulation 
of mass contained in the streams. In addition, after ~ 350 time units, there are slower 
oscillations (time scale ~ 20-30 time units) of the enclosed mass caused by the eccentrically 
orbiting density lump. 

In Figure 3, we display the time-averaged surface density and accretion rate over two 
intervals: t = 250-350 (ATi) and t = 350-450 (AT2). We find there is only minor change in 
the averaged surface density and accretion rate of the disk region < 2.5a over ATi and AT 2 . 
If we define the unit of surface density S = p a, the peak density grows in time, but only 
slowly, rising from ~ O.6X0 to ~ 0.7So- Its position r p also changes, but similarly slowly, 
gradually moving from ~ 2.5a to ~ 3a. The shape of the surface density profile likewise 
hardly changes in time. On the inside, the disk is truncated exponentially, (S(r < 2a))t ~ 
E exp(3.8r/a — 8.4) for both intervals, where ( ) t represents an average over time. Toward 
larger radii, the surface density is oc r~ 2 , the predicted profile of a 'decretion disk' (Pringle 
1991). All that changes is that the region where S(r) oc r~ 2 extends farther outward at later 
times. 
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Averaged over the quasi-steady period, the accretion rate at the inner boundary is 
M ~ 0.018(GMa) 1//2 £ - However, the disk evolves throughout the simulation at r > 3a. We 
see gas inflow rates as large as ~ 0.06(GMa) 1 / 2 S outside the region ~ 3a during both time 
intervals. As a result, the surface density continuously build up around this radius. 

Many properties of these disks can be expected to scale with the mass near the surface 
density peak. It is therefore convenient to define a 'disk mass' = 7rr 2 (E p ) t , where r p = 3a 
and (£ p )t = 0.65So- In terms of this mass (and the natural unit of time for the simulation), 
the mean accretion rate through the inner boundary can be rewritten as 

(M) = 1.0 x l(r 3 M d ft bill . (11) 

The surface density of such disks is very uncertain. A sense of scale, however, can be gleaned 
from translating eqn. 11 into the luminosity that would be produced if the accretion were 
converted into radiation at customary black hole efficiency (10%); in Eddington units, it is 
L/L E = 0.036Mg 7 ' 

a 0A T p^ where we have scaled to a total black hole mass of 1O 8 M , a 
binary separation of 0.1 pc, and a disk column density whose Thomson optical depth is t p . 

3.2. Angular Momentum Budget 

Unlike a steady accretion disk in a point-mass potential, where the internal stress simply 
transports angular momentum outward and thereby drives an inflow, the binary consistently 
interacts gravitationally with the circumbinary disk by torquing the surrounding gas. The 
angular momentum delivered through these torques is transported by two mechanisms: MHD 
stresses, mostly due to MRI-driven turbulence; and Reynolds stresses associated with coher- 
ent gas motions in the gap region. We have checked that the code's numerical shear viscosity 
is negligible compared to the transport mechanisms discussed here. In this section, we first 
investigate the mechanisms of angular momentum transport and the properties of the binary 
torque, and then discuss the balance between binary torques and the stresses. 

3.2.1. Maxwell and Reynolds Stresses 

In the first panel of Figure 4, we show the vertically-integrated and time- and azimuthally- 
averaged stresses and their sum as a function of radius. Within r < 3a, the character of 
these stresses is very similar in our two averaging periods, A7\ and AT2. The Reynolds 
stress shows big fluctuations radially. Its first peak, at r ~ 1.7a, is ~ 0.35Soc 2 . It dips at 
around 3a to nearly zero and then rises up again slowly. The ratio of the first peak to the 
second peak is ~ 10. However, the Maxwell stress has a flatter radial profile. It is roughly 
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~ 0.1S Cg on average and varies only slowly over a range of at most a factor of two. It 
diminishes within the gap and decreases smoothly towards larger distance. Comparing the 
two, we find the Reynolds stress exceeds the Maxwell stress by a factor ~ 4 at r ~ 1— 2a. 
On the other hand, beyond the gap region, the MHD stress always dominates the Reynolds 
stress by a factor > 3. Within one binary separation, the Maxwell stress also exceeds the 
Reynolds stress, mainly because the Reynolds stress falls with the lower gas density near the 
boundary. The total internal stress therefore follows the Reynolds stress at r < 2a, while 
it is close to the magnetic stress outside that region. Similar results have been reported in 
previous studies on gap formation in protoplanetary disks with an embedded planet using 
MHD simulations (Nelson & Papaloizou 2003; Winters et al. 2003). It is not a surprising 
outcome as the disk with a planet is just a special case of a circumbinary with an extreme 
mass ratio. 

In the language of the «ss parameter of Shakura & Sunyaev (1973), we can also measure 
the stresses in units of the pressure. In the right-hand panel of Figure 4, we present the time 
averaged stresses in these units. The most significant distinction compared with the absolute 
stresses is that the stress ratios increase steeply inward due to the low surface density in the 
gap. The total stress ratio peaks at ass — 11 at r = 1.1a. Its value in the disk body is 
almost two orders of magnitude smaller. 

Our simulation shows that the circumbinary disk possesses highly time-dependent struc- 
tures in the horizontal plane, making it important to look at the instantaneous spatial distri- 
bution of the stresses in addition to the time-averaged values. Vertically-integrated snapshots 
of both Reynolds stress (top left) and Maxwell stress (top right) taken at t = 305 time units 
are shown in Figure 5. On top of each plot, we also draw the surface density with 15 contour 
lines, logarithmically spaced from the density floor value to the density maximum. Nearly 
all the Reynolds stress is confined within r < 3a. The most interesting point is that we find 
relatively large stress in a single gas stream extending from ~ l-2a (see the red stripe in 
the top left panel), where the gas gains angular momentum (Sv^ > 0) and is being pushed 
out by the binary. The gas being kicked thus goes on an eccentric orbit, creating negative 
stress at r ~ 2-3a and — 7r/4 < <fi < n/2 (6v r > and 8v^ < 0). In the same range of 
radii but a different azimuthal location (tt/2 < < 57r/4), the gas falls back toward the 
binary, creating positive stress (both Sv r and Sv^ < 0). This bi-symmetric property provides 
a near cancellation so that only the stress within the stream contributes significantly to the 
radial profile of the Reynolds stress. On the other hand, the Maxwell stress is more evenly 
distributed in the disk body, although the typical magnitude is a factor of 4 less than the 
typical magnitude of the Reynolds stress in the inner disk. In the gap, the Maxwell stress 
is strongly positive in the streams, where the field lines are collimated by the gas flow, but 
negative elsewhere. 



-14- 



To illustrate how the stresses depend on the disk height, we present azimuthally-averaged 
snapshots at the same time in Figure 5. We find the Reynolds stress is well concentrated 
in the midplane. Its strength decreases to about half the midplane value at only one scale 
height away from the midplane, and about one tenth at two scale heights. However, the 
Maxwell stress is less sensitive to the altitude. It appears to be mainly confined within four 
scale heights, but also shows positive contributions at higher altitude near the pole due to 
field buoyancy. 



Having discussed the mechanism for angular momentum transport, let us now consider 
the torque exerted on the disk by the binary. We plot the time-averaged torque density and 
the integrated torque in Figure 6, where the torque is calculated approximately using the 
surface density instead of the local density itself. The definition follows MM08: 



where $ = $(r, <f>) is the binary potential in the midplane. 

The first panel of Figure 6 shows the averages of the local torques for the two intervals 
ATi and AT2 are very similar, which again indicates that the inner part of the disk reaches a 
quasi-steady state. The local torque density is positive at r ~ l-2a and peaks at r ~ 1.5-1. 6a 
with f? ~ 0.03-0.04GM£ as the peak value. It then goes negative at r ~ 2-2. 4a, reaching 
its first negative maximum at r ~ 2. 0-2. 2a, but with a smaller magnitude ~ — 0.02GME . 
Similar to what was previously found in MM08 and Cuadra et al. (2009), we find the torque 
density oscillates around zero but damps quickly toward larger distance. In addition, the 
torque density for r < a is negative because the gas within that region is advanced in phase 
(greater angular frequency) with respect to the binary, an effect that would appear in purely 
hydrodynamic treatments as well as MHD. 

The integrated binary torque is displayed in the top right panel of Figure 6. The 
total binary torque exerted on the disk is T(oo) — 0. 011-0. 013GMaSo. In order to make 
a comparison with the torque found in MM08, some normalization is needed. First, we 
normalize the torque to the surface density at the radius where the local binary torque 



3.2.2. Binary Torque 




(12) 



is the local torque, and the integrated torque is 




(13) 
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peaks. We denote that radius as r torq . In our MHD simulation, the first positive peak 
of the torque takes place at r torq ~ 1.6a while it is ~ 1.7a in MM08. In terms of this 
normalization, T(oo)/E(rt or q) ~ 0.07GMa in the hydrodynamic simulation of MM08. In 
our MHD simulation, we find that it is ~ 0.12GMa (averaged in the quasi-steady period), a 
factor of 1.7 greater than in the hydrodynamic and a- viscosity treatment. 

However, as the peak of the torque is located in the gap region, where the surface density 
rises rapidly with increasing radius, this means of normalization is unreliable. We therefore 
turn to a different method, normalizing the binary torque to the disk mass. With this 
method, the normalized torque is T(oo)/Md — 6.5 x 10~ A GMa~ 1 in our MHD simulation. 
However it is only ~ 4.5 x lO~ 5 GMa~ 1 in MM08, ~ 14 times smaller. Just as for the mass 
accretion rate, there is a natural set of units for the total torque: it is convenient to describe 
it in terms of how rapidly the torque removes the binary's orbital angular momentum. In 
that language, 

T(oo) = 2.6 x 10- 3 j bin M d ft bin , (14) 
where j b ; n — (GMa) 1 / 2 /A is the specific angular momentum of the binary. 

Alternatively, as r p is roughly around < 3a for both MHD and hydrodynamic sim- 
ulations, we can write the binary torque as some factor times GMa£ p , and then we have 
T(oo) ~ 1.9 x KT 2 GMa£ p in our present simulation, while MM08 gives ~ 1.4 x 10~ 3 GMa£ p , 
~ 14 times smaller. The largest part of this contrast can be explained by the fact that in 
the MHD treatment S(r t0 rq)/S p ~ 0.15, a factor of ~ 8 greater than that of MM08. 

This contrast can also be understood in terms of the angular momentum budget. The 
gravitational torque per unit mass approximately balances the torque per unit mass due to 
internal stress near the disk inner edge. Therefore, the binary torque should increase linearly 
with the effective a. In our case, a ~ «ss ^ 0.2 in the disk body at r ~ 2-5a, a factor of 
> 20 greater than the constant a = 0.01 assumed by MM08. 

In order to show how the local binary torque evolves with time, we also plot in Figure 6 
the specific torque density (bottom left), i.e., the absolute value of the ratio of torque density 
to surface density. By dividing by the surface density, we remove the effects due to redistri- 
bution of the surface density. We find the specific torque slightly shifts to larger radius and 
diminishes its amplitude from A7\ to AT2. In the bottom right panel, we plot the history 
of the total torque smoothed over a short time span ~ 2T bin , where T bin = 27r/f2 bin is the bi- 
nary period. Before smoothing, the total torque fluctuates very rapidly (~ 1.5f2 bin ) between 
~ — 0.05GMaS and 0.06GMaS in the quasi-steady period. The power spectrum of this 
fluctuation is dominated by a peak at ~ 1.5f2 b i n , approximately the beat frequency between 
the orbital frequency of the matter being torqued most strongly (at the peak of dT/dr, i.e. 
r ~ 1.6 -1.7a) and 2f2 bin . After smoothing, we find the binary torque has a sudden rise 
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around t = 60 which is due to the initial disk filling the small radius region. The torque 
then levels out between t = 100-200. Once in the quasi-steady stage, the torque shows a 
slowly decreasing trend, < 30% in fractional terms over the final ~ 300 time units. It also 
appears to oscillate at the disk orbital frequency (r ~ 2-3a) during the late time, a fluctu- 
ation caused by the eccentric movement of the density lump. As both the disk eccentricity 
and the density of the lump grow with time, the torque normalized by the time-dependent 
EpT-p simply reflects the increasing E p and r p . 



3.2.3. Torque Balance 

For a steady circumbinary disk, the angular momentum at a given radius should not 
fluctuate too much over time. Therefore the angular momentum deposition from the binary 
must either be transported outward by Maxwell and Reynolds stresses or advected towards 
the hole. In order to test this idea, we perform an estimate of the differential angular 
momentum flux averaged over a short time interval between t = 300-320 using the high 
time-resolution rerun. Based on the conservation of angular momentum, we write the angular 
momentum budget for a ring of disk between (r — Ar, r + Ar) as 

— + [(r + Ar) 2 F(r + Ar) - (r - Ar) 2 F(r - Ar)] 

7T dT 
= — [(r- Ar) 2 W^(r- Ar) - (r + ArfW r(f) (r + Ar)] + — , (15) 

where J(r, t) = 27rr{pv ( f ) r)L z is the shell-integrated angular momentum, F(r, t) = {pv r ){pv ( f ) )L z / (p) 
is the mean angular momentum flux due to mass motion, and W rt j,{r, t) = + R r( f, is the 
total internal stress. We write the conservation law this way so that a negative differential 
flux due to advection (the second term on the left hand side of equation (15)) indicates a net 
angular momentum inflow, while a negative flux due to internal stresses (the first term on 
the right hand side) means angular momentum is removed and transported to larger radii. 
We calculate the time averages of all these terms, and plot them in Figure 7 as functions of 
radius. In its y-legend, we use dFj/dr to represent all the terms in the equation, which are 
all radial derivatives of one or another variety of shell-integrated angular momentum flux 
(Fj). 

During this time interval, the local angular momentum stays steady, with only very 
small temporal variations from the averaged dJ/dt (solid black curve), as expected for a 
quasi-steady accretion disk. The binary torque (blue dashed curve) at r ~ l-3a is almost 
balanced by the differential flux due to the sum of Reynolds stress and Maxwell stress (cyan 
dashed curve). In other words, the circumbinary disk within that region reaches a balance 
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between binary torque and internal stresses. As a result, the angular momentum flux carried 
by accretion flows oscillates about zero (magenta curves within l-3a). Clearly, the internal 
torque due to the Reynolds stress (green dash-dotted curve) dominates the Maxwell stress 
(red dash-dotted curve) at r ~ l-3a, which suggests that most of the angular momentum 
dumped by the binary at r ~ 1-1. 8a is transported outward via the Reynolds stress, while 
the angular momentum drawn from the binary at r ~ 1.8-2. 4a is compensated by the torque 
of the Reynolds stress. Finally, we comment that torque balance is not reached at r < a 
(where negative binary torque dominates) or r > 3a (where internal stress dominates). By 
summing all the other torques, we find the differential angular momentum flux that must 
be advected inward in those two regions (magenta curve). Directly computing the second 
term on the left hand side of equation (15) yields the (very similar) black dashed curve in 
Figure 7. 

To sum up, we have found three principal conclusions about the angular momentum 
budget. First, the Reynolds stress dominates the Maxwell stress in the gap region, while the 
latter exceeds the former outside that. Secondly, the distribution of both stresses follows 
the gas streams and are vertically confined: within two scale heights for Reynolds stress 
and four scale heights for Maxwell stress. Third, we find the radial profile of the binary 
torque is similar to that of previous hydrodynamic results, but the amplitude is different: 
normalized by the surface density at the location where local torque peaks, the binary torque 
in the present MHD simulation is twice as great as in previous hydrodynamic calculations; 
however when normalized by the disk mass, our torque is an order of magnitude greater 
than that found in MM08. Lastly, we find that in the quasi-steady state, the binary torque 
is roughly balanced by the torque due to the Reynolds and Maxwell stresses in the inner 
disk (r < 3a). Most of the angular momentum dumped by the binary in the gap region is 
transported outward by the torque of the Reynolds stress associated with the gas streams. 



3.3. Disk Eccentricity 

3.3.1. Evolution of the Disk Eccentricity 

Although we begin our simulation with an axisymmetric disk configuration, which means 
there is zero eccentricity at the start, the disk eccentricity grows throughout our simulation 
(e.g., compare the disk surface density snapshots in Figure 1). 
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To quantify the eccentricity, we define the local eccentricity as 

i (16) 

where ( ) is the shell average defined as in equation (3). We show this local eccentricity in 
the space-time diagram in Figure 8 (top panel). The local eccentricity of the gas in the gap 
(r < 2a) rises significantly during the period t ~ 100-200. The increasing e(r) in the gap 
simply reflects the fact that the dynamics there are strongly non-Keplerian. The eccentricity 
in the main body of the disk (r > 2a) grows rapidly over the period t = 100-300. The 
typical value increases from much less than 0.01 throughout the disk body at t — 100 to 
~ 0.08 between r = 2 -3a at t — 300. It keeps growing slowly after 300 time units, and 
at the end of the simulation a ring of eccentric disk forms at 2a < r < 3.3a with typical 
e(r) ~ 0.1. A more detailed picture of time- variation in the eccentricity is shown in the 
bottom panel of Figure 8, where we show its evolution from t = 300 to 320. For clarity, 
we use the data from the high time resolution simulation, and do not take the temporal 
smoothing as in equation (16). We find strong radial extending structures which connect the 
gap with the disk body. These features bend at ~ 1.5a pointing to later times at both ends, 
which suggests propagation of eccentricity from the gap to both the disk and the binary. 
The pattern repeats once every half binary period, indicating a stream-related origin. We 
will discuss the stream effects on the disk eccentricity in § 4.3.3. 

In Figure 9, we present the radial distribution of the disk eccentricity averaged over 
three time intervals: t = 150-250 (solid black), A7\ (red dotted) and AT2 (blue dashed). 
The distribution curves in the disk body (r > 2a) shift toward larger radius and greater 
eccentricity as the disk evolves in time. Outside the peaks around 2a, the slopes of the 
distribution curves are constant. The radial dependence of e(r) in the disk body can be 
roughly described as oc exp [— 1.3(r/a)]. 

Interestingly, the eccentricity distribution found by MM08 was quite similar to ours in 
both amplitude and shape, despite that calculation's very different internal stresses, initial 
surface density profile, and duration. In our own 2-d hydrodynamical simulations, we have 
likewise found qualitatively similar eccentricity growth. However, the eccentricity profile 
in those simulations matches the MHD simulation (and MM08) only at a specific time 
(t ~ 800fi W n)- From these comparisons, it is clear that the primary mechanism for driving 
eccentricity must be hydrodynamic response to gravitational forcing by the binary, but its 
quantitative development can be affected both by disk physics (e.g., MHD stresses) and initial 
conditions (e.g., the initial surface density profile). We will discuss specific mechanisms at 
greater length in § 4.3.1. 

The space-time diagram suggests we may take the radial average of the local eccentricity 
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for the disk body between r = 2a and r = 4a as a measure of the disk eccentricity. In 
Figure 10, we show the averaged disk eccentricity e<iisk constructed by taking volume averages 
from r = 2a to r = 4a and time averaging over one binary period: 



r.t+T bi: 

Cdisk = tf, — / dt 



J a dr J pv r e l ^r 2 sin 6d6dcp 



bin Jt J 2 ° dr J pv^r 2 sin OdOdcp 



(17) 



At early times (t < 100), Figure 10 shows the transient phase of a non-equilibrium disk in 
a non-Keplerian potential. During this phase, the noise develops significantly, especially in 
the low density region of the initial configuration, and all modes undergo dramatic growth. 
After the transient, the eccentricity of the disk undergoes exponential growth from t = 100 
to t — 250. By fitting the eccentricity curve during that interval, we find the growth rate 
is 7e — 0.018f2bin, where the subscript e denotes eccentricity. The growth rate of disk 
eccentricity slows by a factor of 4 after ~ 250. This suggests some nonlinear mechanisms 
come into play to limit and perhaps eventually saturate eccentricity growth. 



3.3.2. Precession of the Eccentric Disk 

To follow changes in orientation of the eccentricity, we define the longitude of the apoapse 
w in terms of the complex phase angle of the disk eccentricity for 2a < r < 4a. Its 
time evolution is shown in Figure 11. Here the angle w measures the angular location of 
the apoapse with respect to the x-direction in the inertial frame. The disk eccentricity is 
vigorously perturbed by the motion of the gas in the gap during the first one hundred units, 
and the arbitrary orientation in that period is partly because the eccentricity during this 
time interval is very small. The angle w then gradually rotates in a prograde manner. By 
fitting the curve between t = 200-480, we find w ~ 3.2 x 10~ 3 f2bi n - Linear perturbation 
theory predicts the precession rate for particles around the binary potential is w ~ Q — k 
for a cold disk, where Q is the angular frequency and k is the epicyclic frequency of the disk. 
To a first order approximation, Q(r) ss fix 1 + § (" ) 2 h+ q )2 (also see equation (2)), and 

k{t) ps Q K 1 — | (^) 2 -jjfgyi ■ Thus we find the prediction is w ~ 0.19(a/r) 7 / 2 fi b in for q = 1, 
or ~ 4.1 x 10~ 3 £\in evaluated at r ~ 3a. The precession rates measured in both our MHD 
simulation and MM08 are close to this prediction. We therefore suggest that the precession 
is mainly due to the quadrupolar potential. 
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3.3.3. Late Time Lump 



There is another disk characteristic closely connected to the eccentricity: the density 
concentration that appears near the inner disk edge at late time. We call it the 'lump' 
(see Figure 1). The lump orbits around the binary at nearly the same orbital speed as the 
eccentric disk. To quantitatively study the asymmetric feature we introduce the m = 1 mode 
strength component 



which is basically the m — 1 component of the Fourier transform of the surface density. We 
plot the space-time diagram of S m= i in Figure 12. The orbiting lump appears as a zig-zag 
pattern whose peak moves in and out between r ~ 2a and 4a after about 300 time units. 
The pattern period, which is also the orbital period of the eccentric disk, is ~ 30 time units. 
We will discuss the lump in § 4.4 and explain it as due to a combination of disk eccentricity 
and the action of periodic streams impacting upon the edge of the disk. 



As shown in Figure 1, gas is not entirely absent from the gap. Consistent gas flows 
launch from the inner edge of the disk and stream toward the binary. The streams not only 
affect the accretion rate of the circumbinary disk, but also the structure of the inner part of 
the disk. 

In a frame that co-rotates with the binary, we find the streams follow roughly the 
static bi-symmetric potential of the binary and pass through the inner boundary about 0.3- 
0.4 radians beyond the Lagrange L2 and L3 points (see Figure 13). This phase offset is 
consistent with the negative average torque at r < a in § 3.2.2. Early in the simulation, the 
two streams have roughly equal strength, i.e., similar density contrast between the stream 
and the nearby gap region. We show the surface density near the inner edge at t — 122 and 
125 in the first two panels of Figure 13, and draw the velocity field vectors measured in the 
co-moving frame on top of it. The snapshots at t = 122 show the inflowing gas is regulated 
by the binary torque as it leaves the inner disk edge and streams towards the saddle points. 
The sign of the angular speed near the inner boundary alternates from one quadrant to 
the next, a direct signal that the binary torque changes its sign from one quadrant to the 
other. In the co-rotating frame of the binary, the gas in the second (near <p = tt) and fourth 
(near = 0) quadrants is pulled toward the closer component of the binary (see t = 122 
snapshots). The tangential velocity results in a Coriolis force, which tends to move matter 
away from the center. It takes about half the binary period for the gas streams in those two 




(m = 1), 



(18) 



3.4. Streams in the Gap 
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quadrants to be kicked out and strongly interact with the disk edge. As shown in the later 
snapshots at t — 125, the velocity is mostly outward at r ~ 1.5a in those two quadrants, and 
there are strong density enhancements near the regions where the velocity gradient is large. 
There is also an equivalent way to think about the transition between t = 122 and t = 125. 
In the inertial frame, the binary almost always rotates faster than the disk; therefore, the 
second and fourth quadrants have positive binary torque. As a result, gas in these regions 
gains angular momentum from the binary and then drifts away from the center. 

However, the disk inner edge changes its morphology once the eccentricity of the disk 
becomes significant. Instead of a pair of streams with comparable size and strength, we find, 
for instance at t = 301, the stream on the right dominates the left. For an eccentric disk, 
the apocenter side of the disk is far from both members of the binary. Thus, the gas stream 
is strongly torqued, and therefore carries more mass, only when one of the members gets 
close to the pericenter of the disk. After half a binary period, the binary members switch 
their phases (in the inertial frame) and a strong stream forms around the other member. We 
illustrate this process with a time sequence of snapshots at t = 301,302, 303 and 304.1 in the 
corotating frame. The stream on the right at t = 301 is gradually pushed away at t = 302 
and 303 and joins the apocenter of the eccentric disk. Meanwhile, the stream associated 
with the other member of the binary strengthens as that mass approaches the pericenter at 
t = 304.1. 

We also measure the fluid and magnetic effects on the dynamics of the streams. In the 
inner disk and the gap, the force ratio |f g + fL|/|V$| is always < 1/3, where f g = — VP/p 
is the force density of the gas pressure gradient and fi, = (V x B) x VB/47rp is the Lorentz 
force density. Because gas pressure and magnetic forces are smaller than the gravitational 
force, the trajectory of the gas is essentially ballistic. Consequently, we neglect the fluid and 
magnetic effects when calculating the interaction between the streams and the disk inner 
edge in § 4.3.3. 

3.5. Field Structure 

We found in the previous subsection that the Maxwell stress closely follows the gas 
streams in the x-y plane and is well confined near to the midplane. Now we show the field 
structure in the gap region, especially in the stream, has similar features. In Figure 14, we 
draw the vertical averaged magnetic field as vectors on top of the logarithmic scaled surface 
density at t = 300.5. The location and length of the arrows show the position and relative 
strength of the field. Clearly, the field is well collimated by the motion of the gas streams, 
with stronger inward field lines in the lower stream (bottom half plane) and weaker outward 
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field in the upper stream (top half plane). The ordered field in the gap therefore produces 
greater Maxwell stress only within the two streams. 

3.6. Time Dependent Accretion 

We present the accretion rate through the inner boundary as a function of time in the 
top panel of Figure 15. The sampling rate of the time sequence is determined by the output 
rate into the history files, which was once per fi^*. There are high frequency outbursts 
(the narrow spikes in the time sequence) and also a low frequency modulation (the dashed 
curve, which is derived by smoothing the time sequence using a boxcar average with width 
15 time units). These two time variations are the most significant modes in the temporal 
profile of M. By constructing the power spectral density of the accretion rate using Fourier 
transforms (middle panel, Figure 15), we are able to identify them: The higher frequency 
mode is at about twice the binary rotation rate, and it is the same rate at which the streams 
are pulled inward by the binary potential; the lower frequency mode is ~ 0.2f2bi n , which 
is the same as the orbital frequency of the lump during the later stages of the simulation 
(as shown in § 3.3). Our result in Figure 15 is largely consistent with previous findings: a 
dominant frequency associated with the binary orbital frequency (2f2bin if 9 = lj MM08; 
Hayasaki et al. 2008; Cuadra et al. 2009; Rodig et al. 2011; Sesana et al. 2011), a low 
frequency component due to the motion of the dense part of the disk (MM08; Rodig et al. 
2011), and another component (in our case at 1.8f2bi n ) created by a beat between the binary 
orbital frequency and the disk orbital frequency (Rodig et al. 2011). The principal contrast 
between our results and this earlier work is that the dominant frequency is 2f2bm rather than 
flbm because the members of the binary have exactly equal mass. 

To further verify the cause of these two modes, we pick out two times (t = 400 and 441 
as marked in the top panel by black arrows) and plot their two-dimensional local accretion 
rates in the inertial frame in the bottom panel of the same graph. There were accretion 
outbursts at both times, but the former is in the valley of a low frequency oscillation, 
while the latter is at the crest. In the snapshots, the color contours represent the vertically 
integrated accretion rate. The black contour lines show the surface density in a logarithmic 
scale from 10~ 4 to 10~ 5 . These lines pick out the location of the streams, while the white 
contour lines show the surface density between 1.4 to 3.0 in a linear scale; these contours 
delineate the lump region. We find the outbursts of M are indeed coming from the infalling 
streams located in the fourth quadrants of both snapshots. However, the local accretion rate 
in the stream from the left panel is much weaker than that from the right. In the right plot, 
the lump is passing the pericenter region. The binary then pulls a gas stream with relatively 
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higher density than other times. Thus we conclude that the time dependent accretion shows 
both high frequency outbursts corresponding to the periodic inflow from the streams and low 
frequency modulation due to the orbiting lump. The low frequency mode becomes important 
only when the lump forms and maintains itself after about t = 350. 

The time averaged accretion rate at the inner edge, normalized with the peak density, 
is ~ 0.028V GMaSp, a factor ~ 40 greater than found in MM08. Interestingly, during the 
quasi-steady period, the accretion rate in the gap is nearly twice the rate measured at r p , 
where the surface density peaks. This contrast explains the slow mass depletion inside the 
gap shown in Figure 2. One of the most important questions about accretion onto a binary 
is the ratio of the accretion rate accepted by the binary to the rate at which mass is fed into 
the disk from the outside. If the latter rate can be estimated by the maximum accretion 
rate in the disk, we find that accretion through the inner boundary of the simulation, and 
therefore presumably onto the binary, is about 1/3 of the accretion rate supplied. This ratio 
is a factor of two greater than previously found by MM08. However, in evaluating these 
numbers one should be aware that the disk as a whole never reaches equilibrium over the 
course of the simulation. That the time-averaged accretion rate as a function of radius is 
far from constant is a symptom of this fact. Thus, our estimate (and for the same reasons, 
that of MM08 as well) should be taken as no better than an order of magnitude estimate 
of this fraction. We plan a more rigorous approach to this problem in future work. Lastly, 
we note that the effects of the excision on the accretion rate are not explored in this MHD 
simulation. However, our hydrodynamic simulations with different cut-outs suggest that the 
rate hardly changes as long as the cut-out size is smaller than < a. 



4. DISCUSSION 

4.1. Binary Torque and Linear Resonance Theory 

Using the results we have on the disk stresses and the binary torque, we want to further 
elucidate the angular momentum transport mechanisms within a MHD turbulent circumbi- 
nary disk. The binary strongly torques the surrounding gas in the gap as the gas forms 
streams; the outward-moving portions of the streams deliver that angular momentum to the 
inner part of the circumbinary disk. By contrast, the torque directly exerted on the disk 
body is much weaker and diminishes rapidly toward larger distance. In that sense, we can 
say that the binary mainly torques the gas in the low density gap; moreover, the magnitude 
of that torque is proportional to the mass of gas in the gap. This fact suggests that the 
binary torque is very sensitive to the basic physics coupling the streams to the disk such as 
the stresses and the equation of state. For instance, scaled with the peak surface density, 
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our MHD turbulent disk provides a factor of 14 times stronger binary torque than that of 
MM08. 

This fact also leads to a strong contrast with the commonly-used linear resonance theory 
of interaction between disks and planets or other orbiting satelites (Goldreich & Tremaine 
1982). To illustrate that contrast, we calculate both the torque density and total torque 
as functions of radius, assuming most of the torque comes from the f2bin : = 3 : 2 outer 
Lindblad resonance at r2/a = (3/2) 2 / 3 ~ 1.3 (e.g., Meyer-Vernet & Sicardy 1987). The 
estimate is shown in the top panel of Figure 6, with the red lines representing 1/4 of the 
predicted values. The linear theory prediction must be reduced by roughly that amount 
to match the measured total torque. In addition to having a larger amplitude than the 
simulation found, the linear theory also predicts a much shorter wavelength for its radial 
oscillations. 

Three distinctions between the circumstances of the MHD circumbinary disk we simulate 
and the assumptions made in linear theory may explain these contrasts. The first is the sur- 
face density profile. Standard linear theory assumes that any surface density gradients have 
length scales much longer than the wavelength of the density wave excited by the resonance. 
However, the surface density profile in the gap is steep enough that (d In E/dr) -1 < 0.3a, 
which is comparable to the first wavelength 5[(m + l)/(3m 2 )] 1 / 3 (.£/7r) 2//3 a ~ 0.5a, where H/r 
is evaluated at the resonance radius, and m = 2 is used here (see Figure 1 and the pressure- 
dominated case in Table 1 of Meyer-Vernet & Sicardy (1987)). Because the surface density 
at T2 is so much smaller than at r ~ 2-3a, one might expect the maximum torque to move 
away from the resonance point toward the inner disk. Secondly, linear theory assumes nearly 
circular orbits because it also assumes q <C 1. However, in our simulation q = 1, and, largely 
for this reason, the orbits of gas within the inner edge of the disk body are significantly 
non-Keplerian: they form streams. Because the torque depends so strongly on whatever gas 
mass is in the gap, this distinction is qualitatively important. Lastly, as the mass ratio of 
the binary of our simulation is unity, the forcing term, which is roughly controlled by the 
ratio of the perturbing gravitational potential to the sound speed, becomes large enough to 
drive nonlinear density disturbances. The longer wavelength oscillation found in our torque 
density might be partly explained as a result of these nonlinear effects (Yuan & Cassen 
1994). 

4.2. Binary Orbital Evolution 

The time-averaged rate at which the binary loses angular momentum to the disk is 
T(oo) ~ 0.012GMaS . Meanwhile, the time-averaged angular momentum flux across the 
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inner boundary due to accretion is Mj m — 0.017GMaS , where j in ~ 0.94V GMa is the 
averaged specific angular momentum at r in . On average, therefore, the binary actually gains 
angular momentum in net. 

However, even though the binary gains angular momentum, it does not automatically 
follow that the binary separation grows. The evolution of the binary also depends on its 
growth in mass. If the orbit remains circular, the shrinkage rate for an equal mass binary is 

M-4 < i9 > 

a J M 

where J = Mjbin = 0.25MV GMa denotes the angular momentum of the binary, and its 
rate of change is J = Mj- m — T(oo). 

After separating the accretion term of J and combining it with the second term in 
equation (19), we have 

h -2 7 ^ + 2^fl-^V (20) 



a J J \ 2 j 

If jin/jbin > 3/2 (here it is ~ 3.75), the shrinkage rate is controlled by the competition 
between the first and second terms. Plugging in our numbers, we find the two terms are 
surprisingly comparable in mag nitude: -0.096^Ga/MS and 0.082^/Ga/ME respectively. 
The net result is orbital compression, but at a rate much smaller than either of the two terms, 
a/a ~ -0.014^/Gq/MEq. 

However, it is also possible that interaction with the disk may induce eccentricity. If so, 
an extra term on the right hand side of equation (19) should be included, making it 

4 24-s£+(^V)^, (21) 



a J M \l-e bi Je b 



Din 



where ebin denotes the binary's orbital eccentricity. In addition, the angular momentum of 
the binary should also be adjusted to account for the eccentricity: J = Mjbin\A — e bm- ^ n 
Appendix A, we show that the rate of eccentricity growth depends on how the accretion 
flow interacts with the interior disk of each individual binary member, an interaction our 
simulation did not treat. However, we also show that the best estimate we can make from 
the data we do have is that there is little evidence for any significant ebin/ebin- Appendix A 
also briefly discusses the possibility of unstable eccentricity growth raised by linear theory; 
unfortunately, we can say little about it because our simulation data do not bear on it, and 
linear theory may not apply in this regime. If, as we did in § 3.2.2, we define a characteristic 
disk mass Md = 7rS p (3a) 2 , and further assume that the binary orbit remains circular, then 
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a/a ~ -8.0 x l(T 4 (M d /M)fi bin . For comparison, MM08 found a rate ~ -0.003a/ 'Ga/ME , 
which can be scaled used the definition of M d to ~ — 3 x 10~ 4 (M d /M)f2 bin , somewhat slower 
than ours. In other words, although we find a torque an order of magnitude larger than that 
of an a-viscosity hydrodynamic model, we also find an accretion rate so much larger that 
the net effect on the binary orbit is reduced to an increase by only a factor of 2.7. 

Another way to understand this result is to note that our T{pd)/M ~ 0.67(GMa) 1 ' 2 , 
a factor of 3 smaller than in MM08. That is, the greater internal stresses produced by 
MHD effects, especially in the gap region, both introduce more mass to that region (thereby 
magnifying the torque) and lead to even greater accretion, which returns angular momentum 
to the binary. Unfortunately, because the cancellation is so close, our calculation cannot be 
definitive in regard to the shrinkage rate of comparable-mass binaries in Nature. Further 
uncertainty comes from the fact that we cannot track the gas flowing into the cutoff region 
with the present model, so we do not know exactly how much energy and angular momentum 
accretion might bring to the binary, yet our sizable accretion rate makes these contributions 
significant. The net outcome in any particular case may therefore depend on the specific 
details (gas equation of state, binary mass ratio, etc.), as all of these can influence j in , M, 
T(oo), as well as the detailed mechanics of how the stream joins the interior disks associated 
with the members of the binary. 



4.3. Disk Eccentricity Growth 

4-3.1. Eccentricity Distribution 

The distribution of eccentricity in radius plays an important role in the evolution of the 
disk. For the eccentricity to grow substantially over time in a large circumbinary disk, it needs 
to be confined in radius. Otherwise, the energy and angular momentum (more precisely, 
angular momentum deficit) associated with an eccentric disturbance are spread over a large 
radius, resulting in a small eccentricity in any one location. We consider here a linear model 
for the eccentricity distribution and compare it with the results of the simulations. 

The linear equation governing the eccentricity distribution and evolution is given by 
Goodchild & Ogilvie (2006, hereafter, GO06). This also allows generation and damping of 
the disk eccentricity. The basic equation and boundary conditions are given in Appendix B. 
For our circumbinary disk, we adopt input parameters taken from the simulations while the 
eccentricity was small and growing exponentially, i.e., the linear regime. We define the late 
stage of the exponential growth phase as the period t = 200 — 250 time units (see Figure 10). 
Surface density S(r) is taken directly from the simulation by averaging over this time interval. 
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The disk inner radius is taken to be r ; = 2a and the outer radius r Q = 10a. We assume the 
same equation of state as in the simulation: isothermal, with sound speed 0.05Qun a - The 
adiabatic index is taken to be 7 = 1 — O.li. The small imaginary part (which introduces slow 
damping) assures mathematical consistency. We model the forces driving the initiation of 
eccentricty as a Gaussian in radius (see equation (Bll)), centered on r c = 2.1a and having 
a width w = 0.05a. These parameters are not well constrained by the simulation. It follows 
that the the results are also insensitive to their values within certain limits. Parameter s c , the 
peak value of the eccentricity injection distribution in equation (Bll), is determined by the 
condition that the growth rate equal that in the simulation, 7 C = —lm(u) ~ 0.018f2bi n , where 
u is the eigenvalue defined in equation (B14). The value is found to be s c = 0.0145af2bin- 

The physically relevant solution on long timescales corresponds to the fastest growing 
mode. This solution was determined by the methods described in Lubow (2010, hereafter 
L10). Figure 16 plots the eccentricity distribution for this mode, along with the eccentricity 
distribution obtained from the simulation. The middle dashed line corresponds to the eccen- 
tricity distribution based on linear theory with the set of parameters described above. The 
distribution matches well the one obtained from the simulation (the solid line). This distri- 
bution closely resembles the distribution for the fundamental free precession mode in the disk 
(the uppermost dashed line). That mode has an eccentricity injection of zero, s c = 0. The 
confinement of eccentricity is due to the effects of differential precession of the binary d r w g 
(see equation (B8) in Appendix B), while pressure effects act to spread the mode. Differen- 
tial precession effectively creates a potential well that leads to trapping of the fundamental 
mode. Further localization is due to the injection of eccentricity, which is concentrated near 
the disk inner edge (solid curve on lower left). The faster eccentricity is injected into a small 
region, the more the eccentricity distribution narrows. For somewhat higher eccentricity 
injection rates, the distribution is further confined (lowest dashed curve). Therefore, we see 
that the eccentricity distribution is a consequence of the excitation of the fundamental free 
eccentric mode in the circumbinary disk and that the eccentricity is confined within a radial 
width ~ a. 

The excitation of the fundamental mode occurs because that mode overlaps well with 
the eccentricity injection distribution. As discussed above, the location and width of the 
eccentricity injection are not well determined by the simulations. On the other hand, the 
linear equation has solutions that match the simulation results only if r c < 2.6a, which 
indicates that the source of the eccentricity should be located near the inner edge of our 
MHD disk. For larger values of r c , there is no value of s c that could provide the growth rate 
found in the simulation for an eccentricity distribution that resembles the simulated one. 
This radius limit has a weak dependence on the width of the injected eccentricity w. If the 
width is doubled to 0.1a, then the limiting radius increases by about 1%. 
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The outer regions of the simulation exhibit an exponential falloff of eccentricity with 
an e-folding length of about 0.8a (see Figure 9 in § 3.3.1). In Appendix C, we derive an 
asymptotic form for the eccentricity distribution (valid for r ^> a) using linear theory. It is 
dominated by an exponential falloff in r with an e-folding length also of about 0.8a. 

To sum up the analysis in this subsection, we find the linear model can explain the 
eccentricity distribution of our simulation. It shows the confinement of the disk eccentricity 
is mainly due to the disk precession and eccentricity growth. The linear model also puts 
constraints on the location of the source for eccentricity growth. Finally, it explains the 
exponential falloff of eccentricity profile. 



4-3.2. Eccentricity Growth Through Tidal Instability 

One possible mechanism for eccentricity growth is a tidal instability. This instability is 
believed to explain the superhump phenomenon in cataclysmic binaries (see Osaki 1996). 
The instability involves the action of the 3:1 eccentric Lindblad resonance in a disk that 
surrounds a star (Lubow 1991, hereafter L91). As pointed out in L91, this instability can 
also be at work in circumbinary disks at radial locations that satisfy the condition 

n(r) = (22) 

where m is an integer that corresponds to the azimuthal wavenumber of the tidal component 
of the potential involved in the instability. In the case analyzed here of an equal mass binary, 
m = 1 is excluded because that component of the tidal potential is absent. For m = 2, the 
instability is associated with the f2bm : = 2 : 1 resonance that occurs at r ~ 1.6a. In 
this model, the m = 2 tidal field interacts with the m = 1 eccentric motions of the disk to 
produce disturbances of the form exp (3i(f> — 2iQ hin t). These disturbances launch waves of 
that form at the 2 : 1 resonance. The waves in turn interact with the tidal field to produce 
a stress that increases the disk eccentricity exponentially in time. 

However, the 2 : 1 resonance lies at r ~ 1.6a, inside the circumbinary disk gap. Further- 
more, the gas motions are far from circular there. Even if the gas followed circular orbits in 
this region, the model of L91 predicts in this case an eccentricity growth rate ~ aQb in /w e , 
where w e is the width of the eccentricity distribution and is ~ a. This very rapid growth 
rate is a consequence of the powerful m = 2 tidal potential for an equal mass binary, the 
dominant component of the tidal field. However, this rate is more than an order of magni- 
tude faster than the growth rate measured in the simulation. It may be possible that the 
effects of the resonance are still felt, but at a reduced level, near the disk inner edge as a 



-29- 



consequence of finite width of the resonance. Thus, tidal instability may play a role in the 
eccentricity growth, but the evidence suggests it is at most a limited role. 

4-3.3. Eccentricity Growth through Stream Impact 

As discussed in § 3.3.1, the space-time diagram reveals a special feature of the growth: 
the propagation of eccentricity from radii close to ~ 2a toward larger radius. In this section 
we examine the possibility that the impact of streams in the gap striking the inner edge of 
the disk is the mechanism of eccentricity injection into the disk. 

During the exponential growth phase of eccentricity, a pair of gas streams are pulled 
in from the disk inner edge by the binary, and then flung back to the disk after half binary 
period. It is easy to show that if the disk is on a circular orbit and free of eccentricity, the two 
streams are not capable of breaking the bisymmetry that is intrinsic for an equal mass binary 
potential. We speculate that a small eccentricity breaks the symmetry by inducing unequal 
strength stream pairs (in terms of both the density and the velocity) and/or asynchronous 
stream-disk interaction. Once the symmetry of stream impact is broken, stream impact 
could amplify the initially small disk eccentricity, in turn increasing the asymmetry of the 
streams themselves. This process could then lead to sustained disk eccentricity growth. 

We find two sorts of evidence that support this mechanism of eccentricity growth. The 
first is based on a special set of simulations we performed. In these, we eliminated part or all 
of the returning streams by enlarging the central cut-off. If the removal of outward streams 
halts eccentricity growth, we may take that as evidence in support of a stream impact origin 
for eccentricity excitation. We do not consider cutoff radii large enough to affect the circular 
orbit region of the disk because otherwise, the result might be equally well interpreted in 
terms of a resonance model (§ 4.3.2). Because these are essentially hydrodynamic effects 
and involve no vertical dynamics, these special simulations were performed in 2-d with no 
magnetic field (see § 2.5 and Table 1 for details of the simulations). 

In the control run (B2D.rin=0.8), the disk is truncated around 2a, similar to B3D. 
In the gap region, there is a pair of streams at t = 300 — 700 while the disk eccentricity 
grows exponentially. A single strong stream appears when the disk eccentricity becomes 
significant. The simulation is terminated at t — 1500, and at that time the disk eccentricity 
is saturated. We find the stream dynamics and the eccentricity growth in the control run are 
very similar to what were found in B3D, demonstrating that hydrodynamic effects alone can 
yield eccentricity qualitatively similar to MHD. Thus, the results based on hydrodynamic 
simulations can provide a clear indication of the role of the streams in eccentricity generation. 
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We display the growth history of the disk eccentricity in Figure 17 (black solid). Here the 
disk eccentricity is calculated using a two dimensional definition of equation (17). Similar 
to the MHD simulation, we find exponential growth from t ~ 300-700. The growth rate in 
that period is ~ 0.017f2bin, very close to that of the MHD simulation (~ 0.018fibin)- We 
then restart from t = 500 with five different values of rj n while keeping all other parameters 
fixed. We terminate the reruns at t — 1500, the stopping time for B2D.ri=0.8. 

In Figure 17, we plot the evolution of disk eccentricities for cases with different n n . 
The case with r in = 1.0a (blue dash-dotted) loses only the tips of the streams, and its 
eccentricity evolution closely follows the control run, with exponential growth continuing 
without any interruptions until it becomes saturated. The disks with larger cut-offs (r in = 
1.3a: green dashed curve; 1.7a: red dotted curve), however, cease exponential growth of 
their eccentricities right after the restart because most of the streams forming in the gap 
are lost in the hole and never get a chance to interact with the disk edge. There are also 
cases in between. When r in = 1.1a (cyan dash dot dot curve) and r in = 1.2a (magenta long 
dash curve), the growth rate of the eccentricity diminishes after the restart. It takes longer 
times for both disks to reach the eccentricity of the control run. We attribute these results 
to partial loss of outward stream motion. 

When ri n = 1.3a, the simulation domain includes the region well inside the gap where 
gas motions are noncircular and are dominated by the streams. In this case, the simulation 
region also extends well inside the 2:1 resonance at r ~ 1.6a whose presence is required for 
the mechanism of § 4.3.2. Consequently, any resonant instability should not be seriously 
affected by this inner boundary change. However, the eccentricity growth rate did change. 
This evidence then favors the stream impact mechanism. 

The other evidence favoring a stream impact origin of the eccentricity comes from di- 
rectly measuring the rate of change of eccentricity due to the stream-disk interaction in the 
MHD simulation. While the eccentricity grows exponentially, any asymmetry in stream im- 
pacts is necessarily small, so the effects of the gas stream are mild and difficult to directly 
measure. By the time t = 300, when the eccentricity growth is slow, the inner portion of the 
disk forms an eccentric ring between r ~ 2-3a with e ~ 0.08. Its pericenter is at r ~ 2a and 
~ 37r/2, and the apocenter is at ~ 3a and <fi ~ tt/2 in inertial frame. The eccentric ring 
slowly precesses around the binary. In this state there is only a single dominant stream. In 
this configuration, the binary interacts strongly only with the periastron side of the disk. At 
this stage, the effects of the stream are strong enough for us to reliably measure its effects 
on the disk. 

The Gauss equations of celestial mechanics (Brouwer & Clemence 1961) provide a way 
to relate the perturbing effects of the stream to the evolution of disk eccentricity. This 
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method has been previously used to calculate the effects of an external perturbing force on 
a ring in the context of dwarf-nova accretion disks (Lubow 1994). The evolution of the disk 
eccentricity is described by 



de yi — e 2 
dt a^n 



O 4cos/ + 3e + ecos(2/) 

Rsmf + S — — 77 

2(1 + ecos j) 



(23) 



where ad is the semi-major axis of the eccentric ring of disk and n = a/ GMj c? A is the mean 
motion of the disk. Quantity / is the true anomaly where the stream impact occurs. R and 
S are the radial and angular components of the disturbance force density, and are defined 

as 

( g) = m.(y-u,) (24) 

m e 

where m s denotes the mass injection rate of the stream, m e is the total mass of the eccentric 
ring, v s and u e are the stream and disk velocities near the impact, respectively. 

Since the impact takes place over a short time, we must use the higher time resolution 
simulation (t = 300-322) in order to estimate the disk perturbation. We carry out our 
analysis using the vertically-integrated two-dimensional data. For the Gauss equation to 
apply, the properties of the eccentric ring, such as e, ad, n, and the longitude of pericenter 
w, need to be nearly constant over the orbital period of the ring. This assumption holds well 
because the timescale for the ring properties to change is much longer. In addition, the mass 
transferred to the ring by the stream impact over an orbital period is considerably smaller 
than the mass of the eccentric ring. We use time-averaged values for those parameters. The 
values are: e = 0.08, ad = 2.5a and vo = 1.48 radians. The stream-disk impact is localized 
in space. We analyze the stream-disk impact over an arc that is defined by r = 2a and 
4> G [— 7r/4, 7r/4]. The quantities / and v s are directly measured on the arc by taking the 
density weighted averages at each time. The mass injection rate is calculated by integrating 
the mass flux along the arc. The instantaneous velocity of the eccentric edge u e is obtained 
from the unaffected matter at a slightly greater radius, which we take to be an arc at r = 2.1a. 
In fact, the disk velocity is not very sensitive to the location we choose, so long as it is taken 
within the eccentric ring. 

Using these parameters, we calculate the rate of change of the eccentricity from equa- 
tion (23). The result is presented in Figure 18. We find that each impact lasts for less than 
one time unit. The impacts induce peak values of de/dt that range between 4.0 x 10~ 3 f2bi n 
and 0.013fibin at each half binary period. Most of the time, de/dt remains close to zero 
because stream-disk impacts are so brief. Taking time averages of the curve, we obtain an 
estimate of the impact-induced rate of change of eccentricity {de/dt) t — 1.5 x 10 _3 fibin- The 
corresponding growth rate 7 e ~ 0.019f2bin is consistent with the growth rate at earlier times 
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in the simulation (~ 0.018fibin)- The growth rate contribution predicted by the stream im- 
pact at this later time may not be quite the same as the impact contribution at the earlier 
stage. Nonetheless, the approximate correspondence is encouraging. Indeed, one question 
that arises is why, given the estimate just made, the growth rate is smaller at this later time. 
We can only speculate that nonlinear effects, significant at this time but not earlier, may 
limit eccentricity growth. 

4.4. Interpretation of the Lump 

The disk after t = 400 shows a significant asymmetry which we call the "density lump" . 
We believe the lump is due to the combined effect of stream impact and disk eccentricity. As 
discussed in § 3.4, only one stream strongly interacts with the binary at any given time once 
the disk is noticeably eccentric. The stream moves outward after gaining angular momentum 
from the binary. Once it reaches the disk inner edge, the stream is shocked and its fluid 
compressed. Its mass is added to the nearby gas mass, increasing the local density. This 
region of enhanced density then moves around the binary at about the local orbital speed, 
decaying over roughly one orbital period due to pressure forces and shearing. One might not 
therefore expect such a small scale density enhancement to become a large scale, long-lasting 
lump. 

We have, however, identified two mechanisms that sustain the lump and foster its 
growth. First, the lump can grow more concentrated via streams coming from the lump 
itself. We name this mechanism stream reabsorption. When the lump approaches percen- 
ter, the binary peels off a gas stream that is denser than streams drawn from lower density 
regions. This relatively dense stream is then kicked back out by the binary torque. Although 
the azimuthal location reached by the returning stream may not be exactly where its start- 
ing point has arrived at this time, given the relatively large azimuthal extent of the lump, 
the chances that the returning stream strikes somewhere in the lump are relatively good. 
Moreover, the forward shock driven by the returning stream into the lump gas compresses 
it, restoring the density loss incurred by shearing during the time the stream passed near 
the binary. The left panel in Figure 19 shows the moment an episode of stream reabsorption 
taking place. 

The second mechanism makes use of streams drawn from regions of the inner disk other 
than the lump. We call it lump feeding, as this channel of development not only enhances 
the concentration of the lump, but also increases its mass. After t = 400, the density 
enhancement orbits eccentrically, and as it follows its orbit toward apocenter, both v r and 
v,], diminish. Meanwhile, streams traversing the gap continue to strike the disk inner edge, 
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creating a new, smaller density enhancement at orbital phase angles behind the lump. The 
orbital speed of this new density concentration is much greater than that of the main lump 
because it is at rather smaller radius. As a result, the newly-created and smaller lump 
can catch up with the main lump and join it before the lump returns from apocenter and 
accelerates. In the right panel of Figure 19, we show an example of this process, pictured 
at a time (t = 422) shortly before an outgoing stream reaches the lump. Three time units 
later, the small enhancement (the green region close to r = 2a at <fi < tt/2) joins the main 
lump (green and red colored area in the second quadrant). The density- weighted velocity in 
the enhancement is (v r , v^) ~ (0.17, 0.66)f2bi n a (we average the stream gas where the surface 
density is greater than 1.5E ). However, the lump has a slower velocity ~ (0.07, 0.52)f2bi n a 
(averaging over locations where S > 2S ). As a result of both the lump feeding and stream 
reabsorption, the density lump grows both in mass and density contrast. 

Another way to illustrate the connection between streams and the lump is to look at 
the space-time diagram in Figure 12. In the diagram, the stream clearly stands out as the 
radially stretched features between r < 1.5a and r > 2a which propagate away from the 
binary with a pattern speed of about 0.2n bin a. This pattern repeats itself at a rate about 
twice the binary frequency. Meanwhile, the zigzag oscillation at late time (t > 400) between 
radius ~ 2a and 4a shows the movement of the density lump as it orbits around the binary 
on an eccentric orbit with a period of ~ 30 time units. Clearly, when the lump passes 
the pericenter, the stream reabsorption process promotes growth in the density contrast, 
driving an increase of the mode strength at r = 2-2. 5a in the ascending part of the zig-zag. 
The ascending part at even larger radius is able to maintain its strength via lump feeding. 
However, once the lump passes apocenter, the lump feeding is limited as the lump is now 
much further away from the stream in azimuthal angle, and thus we find the descending part 
at r > 2.5a is less affected by the stream features. As the lump approaches pericenter again, 
another cycle of stream reabsoprtion and lump feeding begins. Note that the blue dots (low 
mode strength regions) around 2.5a between the zig-zags actually show the moments that 
the newly kicked out stream reaches the radial coordinate (but is distant in azimuth) of the 
lump, so that it smooths the m — 1 mode at that radius. The zig-zag feature's growth in 
both strength and radial range demonstrates that the density of the lump increases gradually, 
and as the disk becomes more eccentric the semi-major axis of the lump's orbit grows. 
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5. CONCLUSIONS 

5.1. Specific Results 

In this work, we have performed the first three-dimensional MHD simulation of a cir- 
cumbinary disk around an equal- mass binary, which in this case follows a circular orbit. The 
main results are: 

1. The disk exhibits a number of nearly steady features. There is a low-density gap within 
r < 2a , an eccentric inner disk at (2-3) a. At early times, there is a pair of gas streams 
that flow into the gap from disk inner edge. They are nearly steady in the corotating 
frame of the binary. At later times, there is a single dominant stream whose mass flux 
is time varying with a period of half the binary orbital period. Parts of these streams 
are torqued so strongly by the binary that they return and impact on the disk inner 
edge. 

2. Some aspects of the disk evolve secularly. At late time, the disk inner edge develops 
an asymmetric density concentration ('the lump') whose mass, density contrast, and 
orbital eccentricity grow steadily. We find the lump is due to a combination of stream 
impact and disk eccentricity. 

3. The disk eccentricity grows exponentially during the time t = 100-250 (16-40 binary 
orbits) with a growth rate of ~ 0.018f2bin- The growth rate then slows down signifi- 
cantly. By the end of the simulation, the disk body at 2a < r < 3a has reached an 
eccentricity ~ 0.1. Its pericenter precesses slowly due to the quadrupolar component 
of the binary potential. Stream impact largely accounts for the eccentricity growth. 

4. Reynolds stress associated with the streams is the leading transport mechanism in the 
gap region, while Maxwell stress dominates in the disk body. 

5. Although the profile of the binary torque is similar to that of previous hydrodynamic 
simulations, its magnitude depends strongly on the magnitude of internal disk stresses. 
Normalized to the disk mass near the density peak, we find a measured total torque 
T(oo) ~ 14 times greater than found by MM08. Because the torque at any particular 
radius is proportional to the gas mass there, the integrated torque is directly propor- 
tional to the gas mass in the gap region, where the binary torques are strongest. This 
mass, always a small fraction of the disk mass, reaches the gap through the action of 
internal stresses in the accretion flow. Relative to the gas pressure, MHD stresses in 
the disk body are about an order of magnitude larger than the "viscous" stresses esti- 
mated phenomenologically in previous work; in the gap, this ratio increases by another 
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order of magnitude. As a result, the density of matter in the gap is also about an order 
of magnitude larger than found in hydrodynamic calculations, and this contrast leads 
directly to the larger torque. 

6. After time averaging, the accretion rate at the inner boundary is ~ 30% of the peak rate 
in the disk body. Compared with MM08 and normalized to the peak surface density, 
the accretion rate at the inner boundary in this simulation is ~ 40 times greater. This 
contrast, like the contrast in the total torque, can be largely attributed to the stronger 
internal stresses due to MHD effects. The time-dependent accretion rate is strongly 
modulated on both the binary orbital frequency (by the stream) and the inner-disk 
orbital frequency (by the lump). 

7. Previous work on the angular momentum budget of the binary has focussed on the 
torque it exerts on the disk. Because we find a substantially larger accretion rate than 
previous calculations, we also find that the binary's gain in angular momentum due to 
accretion can substantially offset its loss by torque. As a result, the estimated binary 
contraction rate a/a ~ — 8 x 10~ 4 (M d /M)(GM) 1 / 2 a~ 3 / 2 is only slightly larger than the 
rate estimated by MM08. 

5.2. Consequences for Orbital Evolution: Md < M vs. Md > M 

These detailed results have a number of more general implications for the evolution of 
circumbinary disks, particularly in the context of supermassive black hole binaries. It is 
generally believed that stellar dynamical effects become relatively ineffective at driving the 
evolution of this sort of binary when its separation falls much below ~ 1 pc (Begelman et 
al. 1980). Recently, there have been several attempts to solve this 'final-parsec problem' 
with stellar dynamics. Although they are potential solutions to this problem, they require 
either special non-axisymmetric stellar distributions (Berczik et al. 2006; Khan et al. 2011; 
Preto et al. 2011) or extra perturbers such as giant molecular clouds (Perets & Alexander 
2008). Given the uncertainty about whether these candidate mechanisms suffice to solve 
the problem, the prospect that angular momentum loss to a surrounding disk may push the 
binary through this barrier is an attractive one (e.g., Ivanov et al. 1999; Gould & Rix 2000; 
Armitage & Natarajan 2002; Escala et al. 2005; Kocsis & Loeb 2008; Schnittman & Krolik 
2008; Cuadra et al. 2009; Dotti et al. 2009; Corrales et al. 2010; Tanaka et al. 2011; Farris 
et al. 2011). 

We have previously described this process in terms of the orbital shrinkage rate a /a. 
Defining the orbital shrinkage time t s hrink = |a/d| , the time taken by torque alone to change 
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the binary orbital angular momentum t tor que = Mj Un /T(oo), and the disk accretion time 



Using the values found in the simulation makes taccAsimnk — 0.6M d /M, i.e., the ratio of the 
accretion time to the orbital shrinkage time is roughly the same as the ratio of the disk mass 
to the binary mass. Indeed, Lodato et al. (2009) argued that, unless the disk mass were 
at least comparable to the binary's secondary mass, it would be ineffective in driving binary 
evolution. Although it is true that a mass comparable to the secondary's mass must pass 
through the inner region of the disk over an orbital evolution time £ s hrink> it is not necessarily 
true that this mass must be there all at once. 

Suppose, for example, that Ma <C M. Because t acc <C £ s hrmk in this case, if the disk 
mass is put in place once and for all, it would certainly be drained long before the binary 
has significantly evolved. On the other hand, one could also imagine situations in which the 
disk is continuously replenished. In that case, the binary orbit could change substantially 
due to this interaction even though the instantaneous disk mass is always much smaller than 
the binary mass; it's just that this process takes longer. Before determining how much time 
is required, it is important to recall that the close subtraction between the torque term and 
the accretion term may make the numerical value of (t a ccAshrmk)/(Mi/^0 rather sensitive to 
specific assumptions of our simulation, e.g., the binary mass ratio and the disk gas's equation 
of state. If the mass accretion rate were reduced by a factor of two while the torque was 
held fixed, the shrinkage rate would increase by a factor ~ 4; conversely, if it were increased 
relative to a constant torque by a ratio > 1.2, the binary orbit would actually expand over 
time. 

These processes could also be affected by the fact that the accretion rate onto the 
binary is a substantial fraction of the accretion rate in the outer disk. The specific number 
we found (~ 30%) is not terribly well-defined because the lack of inflow equilibrium in the 
outer disk makes the denominator in that ratio very uncertain. However, a reinterpretation 
of that figure as, more vaguely, tens of percent, nonetheless has significant implications. If 
the accretion rate in the outer disk were sufficient to supply an AGN, even the fraction 
leaking through the disk's inner edge would still be large enough to fuel AGN activity, even 
if somewhat weaker. As a result, the inner edge of the disk would be illuminated in much the 
same way AGN are known to illuminate the inner edge of their "obscuring tori" (Antonucci 
1993). Low density gas in the gap region would then be strongly heated by absorption 
of the AGN continuum, and much of it driven outward in a wind (Krolik & Begelman 
1986; Balsara & Krolik 1993). Such a wind would lead to a reduction of the rate at which 



t 



acc — 



M(j/ M, we can rewrite this rate as 




(25) 
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mass is captured by the members of the binary. As a result, the AGN luminosity would 
be smaller, and an equilibrium in which only a fraction of the accretion rate through the 
disk inner edge is captured by one or the other of the black holes might be achieved. An 
immediate consequence would be an increase in the ratio of torque to angular momentum 
captured by the binary, and therefore a factor of several increase in the shrinkage rate. Of 
course, a qualitative change in the equation of state of gas in the gap might well lead to an 
order unity change in the torque, as well. At the same time, radiation forces associated with 
illumination of the disk body can thicken it (Pier & Krolik 1992; Krolik 2007; Shi & Krolik 
2008). It is unclear how such a situation (effective vertical gravity weaker than zfl 2 ) would 
alter saturation of the magneto-rotational instability, but more rapid inflow is one plausible 
consequence. 

The opposite case is also worth considering, in which 3> M, so that t acc 3> tshrink- If 
this is so, then even a one-time deposit of mass in the disk might drive strong evolution in 
the binary orbit. However, when the binary separation decreases, the position of the strong 
torques moves to smaller distance from the center of mass as well, raising the question of 
how much mass might be there to receive those torques. We can estimate that time within 
the disk body of our simulation via tmflow(^) = j r dr'/(v r ) Pt t, where (v r ) Pt t represents the 
time- and density-weighted shell- averaged inflow velocity. In a steady-state disk around a 
point-mass, t acc would match ti n fl ow at r = r p . However, to the degree that binary torques 
retard accretion, t acc exceeds the inflow time measured at the surface density peak, which is 
£inflow( r P ) ~ 5 x lO 2 ^*. At sufficiently large radius, t in Q ow (r) will nearly always be greater 
than t acc because tmflow(^) oc r 3 / 2 (r / H) 2 ; only a rapid outward flare in the disk thickness 
could alter this conclusion. It has long been thought that if t in fi ow (r p ) ~ t acc 3> Shrink, the 
disk and the binary would decouple once the binary shrinks by a factor of a few because too 
little gas would be found close enough to the binary to feel the torques; if orbital compression 
depends solely upon interaction with a circumbinary disk, the shrinkage time could therefore 
never be much shorter than t- m f[ ow (r p ) (e.g., Gould & Rix 2000; Armitage & Natarajan 2002). 

However, the dynamical behavior seen both in our simulation and in others hints that 
this assumption might not be valid. The streams occurring in our simulation travel inward 
on a timescale ~ ^bhu which is <C inflow Although a smaller binary separation would make 
the trajectories of the streams at radii not far inside the disk edge more nearly angular 
momentum-conserving, and therefore make it more difficult for them to penetrate to much 
smaller radii, the internal magnetic stresses within the streams might be strong enough to 
replace the binary torques. Indeed, we found that the nominal "a" in the gap region in 
our simulation was one or two orders of magnitude larger than in the disk body. There is 
also another mechanism, commonly seen in global MHD simulations of disks around point- 
masses, that might contribute in this way. As emphasized by Guan, Hawley & Krolik (2011), 
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it is easy for disks far from inflow equilibrium to drive mass-inflow rates much greater than 
those expected on the basis of inflow equilibrium estimates. For example, when an annular 
disk with a seed poloidal magnetic field is created, orbital shear rapidly creates toroidal field 
from radial field. The —BrB^ stress acting on the low density matter near the disk's inner 
edge rapidly removes its angular momentum, propelling it inward. As the binary separation 
becomes small compared to r p , the gravitational potential comes to resemble a point-mass, 
and similar dynamics might be expected. Because the exertion of a sizable torque requires 
only a small fraction of the disk mass to orbit in the gap region, where the torque per unit 
mass is greatest, a strong torque might continue even though the bulk of the disk mass 
remains far from the binary. 

A massive disk is also likely to be subject to self-gravity, a mechanism not treated in our 
simulation. If M^/M > (c s /v ox b), self-gravitational instability can be triggered (Goldreich 
& Lynden-Bell 1965). The evolution of this instability depends on the ratio of the cooling 
time to the dynamical timescale. If the ratio is greater than unity, the disk can maintain 
marginal stability, and the self-gravity would add extra stress to facilitate the accretion 
(Gammie 2001; Lodato & Rice 2004, 2005). If the cooling time is short, which is likely to 
be the case for disks around massive black holes, the disk is likely to fragment and much 
of its mass may be transformed into stars (Shlosman & Begelman 1989; Gammie 2001; 
Nayakshin 2006), reducing the gas surface density (Lodato et al. 2009). The consequences 
for interaction with the binary remain somewhat unclear. In sufficiently massive disks, stars 
might form with orbits taking them sufficiently close to the binary that the summed torques 
from many individual star-binary encounters could continue to shrink the binary; in effect, 
the stellar loss-cone is repopulated by local star formation (Cuadra et al. 2009). In less 
massive disks, stars might form only in the disk body. If most of the disk mass is converted 
into stars, the gas mass available to fill the gap would be reduced, possibly leading to both 
a smaller torque and a smaller accretion rate. On the other hand, gravitational torques due 
to non-axisymmetric fluctuations in the stellar density might create stresses strong enough 
to overcome the reduction in gaseous disk mass. (Binney & Tremaine 2008). Given all these 
complicated possibilities, it is clear that further study of self-gravitating circumbinary disks 
will be necessary before any strong conclusions can be drawn on the fate of massive disks. 



5.3. Heating From Stream Impact 



The work done by the binary on the streams can heat the inner edge of the disk through 
stream impact. The total energy per unit time delivered by the binary torque can be esti- 
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mated as 

^t orq = ^binT(oo) ~ 6.5 x lO^n^a 2 M d . (26) 
For disks around massive binary black holes, 

L torq ~ 1.9 x 10 38 M 3/2 a ^ /2 iV 24 ergs s~\ (27) 

where Mq = M/1O 6 M is the binary mass, ao.i = a/0.1 pc the separation, and N24 = 
Nn/1.7 x 10 24 cm -2 is the column density of the disk. Thus, the expected luminosity is well 
below a typical AGN luminosity, and would therefore be difficult to detect unless its power 
were emitted in a small number of lines. Similarly, we can get the heating rate for disks 
around stellar binaries: 

L torq ~ 7.9 x 10 33 M 3/2 a- 5/2 Mi,-i ergs s" 1 , (28) 

where Mq = M/M is the mass of the binary in solar masses, aio = a/10 AU is the binary 
separation in units of 10 AU, and M^-i = M&/0.1 M is the disk mass in units of 0.1 M . 
We note that, with our ideal MHD model, this luminosity only sets the upper limit of the 
possible heating rate due to the binary torque. In reality, disks outside 10 AU might possess 
low-ionization dead zones (Gammie 1996) which suppress the MHD turbulence and diminish 
the mean effective ass- Based on the accretion rate observed in disks around young stars, 
which is typically ~ 10~ 8 M yr -1 (e.g., Hartigan et al. 1995; Gullbring et al. 1998) for T 
Tauri stars, we expect the actual luminosity will be two to three orders of magnitude smaller 
than estimated in equation (28). Because stream impact is periodic with the period of the 

]/2 3/2 

binary, the resulting luminosity might be modulated with a period ~ 16M a 1Q yrs. 



5.4. Implications for the "Last Parsec Problem" 

In the black hole binary context, circumbinary disk studies have been in part motivated 
by a search for mechanisms to solve the "last parsec problem", the expected slow-down in 
orbital evolution by stellar encounters when the binary separation becomes < 1 pc. Our re- 
sults have, in one sense, quantitatively weakened the constraints on this proposed solution: 
per unit disk mass, our nominal result is an orbital shrinkage rate a few times faster than 
previously thought. However, in a number of other ways, we have raised potential compli- 
cations. One is that the nominal shrinkage rate is very sensitive to parameters because the 
larger accretion rate per unit disk mass we find leads to a near cancellation between the bi- 
nary's loss of angular momentum through torques on the circumbinary disk and acquisition 
of angular momentum through accretion. Another stems from the likely back-reaction on 
both the disk proper and the matter in the gap caused by AGN illumination as a result of 
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accretion onto the binary. Still another is the question of how the asymmetry in the disk (the 
"lump") affects dynamics driven by disk self-gravity when the disk is sufficiently massive. 
For all these reasons, our work has enriched, rather than settled, the question of whether, or 
to what degree, circumbinary disks alleviate the "last parsec problem". 

5.5. Future Directions 

Finally, we acknowledge that, as the first MHD simulation of a circumbinary disk with 
order-unity mass ratio, this work is limited in several aspects. 

Most of all, it would be interesting to know the long-term behavior of the lump. The 
reason why we ended the simulation when we did was that matter accumulates more rapidly 
than magnetic field in the lump. Consequently, its internal Alfven speed diminished, and 
the ability of the simulation code to resolve adequately the MHD turbulence along with it. 
Such a situation leads to an artificial weakening of the magnetic field. To follow properly 
what actually happens in the lump will require significantly better spatial resolution. 

There are also numerous different parameters of the disks that need to be explored, 
such as the mass ratio and the disk thickness. For example, a mass ratio different from unity 
might result in noticeable changes in both the structure of the disk and the leakage rate 
because the secondary would come closer to the disk inner edge. A smaller thickness would 
likely reduce the inflow rate in the disk body. Such a change might diminish the ratio of gas 
mass in the gap to disk mass, but the considerations discussed above about MHD dynamics 
specific to the gap might counteract this effect. 

We have also assumed that the disk orbits precisely in the binary orbital plane, and in 
a prograde sense. Inclined orbits might well occur, and their interaction with the binary's 
quadrupole potential will lead to precession of the inner disk's orbital plane. Coplanar, but 
retrograde, gas orbits eliminate the resonance structure characterizing the linear theory of 
binary-disk interaction and could lead to other qualitative effects (Nixon et al. 2011). 
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APPENDIX 



A. Binary Eccentricity Evolution 



In this section, we estimate the growth rate of the binary's eccentricity based on a 
simple model of accretion and energy exchange between the disk and the binary. The orbital 
eccentricity e^n can be written in terms of the total orbital energy E, angular momentum 
J, and total mass M: 



where \x = m\mij {m\ + m 2 ) is the binary's reduced mass, m,\ and m-i are the masses of the 
binary components, and M = mi + m 2 . Taking the time derivative of equation (Al), we 
find ew 



Note that if the mass ratio mi/m 2 = 1 (as in our simulation), e^ n is independent of any 
change in mi/m 2 . 

In order to estimate the rate of change of the orbital energy E, we must understand 
how the binary and circumbinary disk exchange energy. The binary can lose orbital energy 
by doing work on the surrounding disk through its time-averaged torque. In addition, the 
accretion flow penetrating the gap ultimately becomes attached to the binary, bringing its 
energy to the binary. However, when the accreting matter begins to orbit around an individ- 
ual member of the binary, we must distinguish how much of its energy becomes associated 
with that motion as opposed to the binary orbit. In addition, there may also be radiation 
losses. Both of these effects could be computed if the binary members were on the grid of 
the simulation, but they were not in the calculation we performed. Consequently, the best 
we can do here is to bound the range of possibilities. 

For the greatest physical insight, it is convenient to group the contributions according 
to whether they are associated with the binary-disk torque or with accretion, i.e., writing the 
terms in the parentheses of equation (A2) as (E/E + 2 J / J) t0 rq + (E/E + 2 J/J) aC c — 5M/M, 
where the subscript 'torq' denotes torque related contributions, and 'acc' the accretion related 
quantities. 

Consider the binary-disk interaction first. When the binary's orbit is exactly circular, 
the work done by the disk on the binary can be written in terms of the azimuthal component 




(Al) 




(A2) 
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of the binary's gravity acting on the disk: 

W = J E(r, <j>) (J^) n hhl rdrd(f> = -T(oc)ft bin . (A3) 

In our simulation, this contribution (in terms of the change in binary orbital energy) is 
~ — 0.012GMaS fibin (also see section 5.3). Because the orbital energy of the binary is 
-GM 2 /(8a), the ratio W/E = 0.096(a 2 S /M)fi bin . The rate at which the orbital angular 
momentum changes due to torque on the disk is — 0.012GMa£ ; hi ratio to the binary 
angular momentum, (l/4)My/GMa, this contribution to J/ J is — 0.048(a 2 £o/M)f2bm- The 
combination (E/E + 2J/J) tor q for these two pieces is therefore zero to within the accuracy 
with which we can do the calculation. 

To describe the energy and angular momentum brought to the binary orbit by accretion, 
we employ the following toy model. Around each of the binary member, the accreting gas 
would presumably form an interior disk: a disk extending out to a distance r& from its 
central mass. We imagine each binary component together with its disk as a 'hard sphere' 
of that radius orbiting around the center of mass of the binary system. When the stream of 
accreting matter reaches the sphere's surface, it sticks to the sphere, adding its momentum 
to the sphere's. It also delivers its potential energy. For simplicity, we will deal only in 
quantities averaged over times much longer than the binary's orbital period, but shorter 
than the timescale of the binary's orbital evolution. 

In this model, the rate of change of kinetic energy in the binary orbit due to accretion 

is 

E k = Mv s ■ v hin , (A4) 

where E^ is the binary's kinetic energy, v s the velocity of the stream, and Wbm the velocity 
of the sphere, i.e., the orbital velocity of the individual binary component. For an equal 
mass binary, the time averaged |"Z?bin| = (1/2) ^/GM/a with respect to the center of mass of 
the binary. From energy conservation, we see that |if s | = A/2e- n + GM/a + GM/r^, where 
e| n is the specific energy of the stream when it hits the 'hard sphere'. This energy could be 
different from e- m , the one measured when the stream crosses the simulation boundary. We 
approximate = 0.3a as the radius of the tidally truncated interior disk for an equal-mass 
binary (Paczyhski 1977). Mass addition from the stream diminishes the binary's potential 
energy by E p = —M(GM/2a). This can be viewed alternatively as the potential energy of 
the mass stream with respect to the other binary member or as twice the potential energy 
per unit mass of the binary as a whole because there is a contribution both from the arriving 
matter and from the deepening of the potential due to the increase in total mass associated 
with accretion. The potential energy of the stream with respect to the binary member it is 
joining is associated entirely with the orbit of the stream around that member and does not 
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contribute to the orbital energy of the binary. Summing the kinetic and potential energy 
brought to the binary by accretion and forming its ratio to the orbital energy, we find 




4^ 
M 



GM/a r d 



(A5) 




where (cos^) is the time averaged angle between v s and Vb in . Assuming e- n ~ e in = 
— 0.81GM/a (its value measured at the simulation inner boundary r = r; n by taking mass 
flux weighted shell averages of the kinetic and potential energy), the accretion contribution 
to the change in binary energy is 

M 

4— (1-I.65(cos0)). (A6) 

Angular momentum is brought to the binary at the fractional rate 

M 

= 3 - 76 m (A7 > 

because the specific angular momentum of the binary is (l/A)yGMa and the specific angular 
momentum of the stream when it crosses the simulation boundary is j- m = 0.94\/GMa. 

The total rate of change of eccentricity (see equation (A2))is then the sum of the accre- 
tion contributions due to energy, angular momentum, and mass acquisition 

!^i^!(6.59<cos*)- 6.52)! (A8) 

We have previously found that M/M = 0.018(a 2 £ /M)ft bin , so that 

6 — = (6.44(cos 8) - 6.36) x 1(T 3 fi bin . (A9) 

e bin 2e 2 V M ) 

In other words, the eccentricity decreases unless {cos 8) is very nearly unity. Because our 
simulation assumes a circular orbit, this result seems surprising because there is no obvious 
reason why (cos 8) cannot be significantly smaller than one. In fact, however, the orbital 
mechanics require (cos 8) ~ 1, as can be seen from the following argument. 

When measuring the specific orbital energy of the stream as it crosses the inner problem 
boundary, we find its azimuthal velocity is ~ 1.21<y/ 'GM/a, almost twice the magnitude of 
its radial velocity (~ —Q.QQ\jGM / a). Thus, the stream's velocity is nearly in the azimuthal 
direction. Moreover, this azimuthal speed is rather greater than the binary orbital speed 
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0.5^GM/a, and, as shown by the six panels of Figure 13, most matter crosses the inner 
boundary at an orbital phase slightly ahead of the nearest member of the binary. Therefore, 
when the stream crosses the boundary, it does so at a small angle and then rapidly catches 
up with the other member of the binary, traveling at almost constant radius. In other words, 
when it strikes that disk, it does so with {cos 8) ~ 1. 

The net result is that the accretion contribution to eun/^bm, as well as the disk interac- 
tion contribution, almost exactly cancels. To within the accuracy of this simulation, these 
effects appear to have at most a weak effect on the binary eccentricity. 

There is, however, a resonant interaction between binaries and their surrounding disks 
that has the potential to drive a linear instability in the binary eccentricity. The resonance in- 
volves forcing by a binary potential component of the form $ m ,z(r, 0, t) = m ,/(r) cos (m0 — fibm^)- 
It is strongest at the 1:3 resonance, where m = 2 and / = 1 (Artymowicz et al. 1991). Using 
the same formalism as previously, the rate of growth of eccentricity can be written as 

dine _ E - 

IT ~ 2^\eY (A10) 

when e C 1. Using E = Q p J with Q p = Q^ n l/m = fibm/2 for the 1:3 resonance, it then 
follows that 

^ ln 6 - (\^^\ 
~dT ~ ~2e 2 /ift bin a 2 ' 

where fi is the reduced mass of the binary. To evaluate this eccentricity growth rate for the 
resonance, we determine the torque J = —T^\ as defined in equation 11 of Artymowicz & 
Lubow (1994). This torque is determined through the use of equations (12)-(14) and (21) 
in that paper. We obtain 

dine 497T 2 m!m 2 a 2 S 

~dT = 16 M 2 If (M2) 
where surface density £ is evaluated near the resonance at r ~ 3 2//3 a and for an equal mass 
binary m\ = = M/2. If the density varies near this radius, then it should be suitably 
averaged over the resonance width of order ~ (H/r) 2 ^ 3 r. There is a caveat regarding this 
result, however: it is derived in the context of linear theory and assumes an axisymmetric 
disk whose fluid follows circular orbits. It is uncertain how much the growth rate might 
change in a highly non-axisymmetric, eccentric disk like the one found in our simulation. 
Moreover, once the binary eccentricity begins to grow, new resonances appear whose effect 
tends to counteract eccentricity growth (Lubow & Artymowicz 1992). Thus, it is hard to 
evaluate the ultimate impact of this possible instability. 



-45 - 



B. Linear Equation for Eccentricity Distribution 

We apply the linear equation for eccentricity evolution given by GO06 that can be 
written as 

id r {h{r)d r E{r,t))+ih{r)E{r,t) + J(r)s(r)E(r,t) 

= J(r)d t E(r,t), (Bl) 

where E(r,t) = e(r,t) exp (izu(r,t)) is the complex eccentricity, for real eccentricity e and 
periapse angle w. 

E is related to the linear perturbations from the axisymmetric circular velocity. For the 
velocity expressed cylindrical coordinates as (u'(r, t) exp (— i<f>), v'(r, t) exp (— i(f>)), we have 
that 

u'{r,t) = irQ{r)E{r,t) (B2) 

and 

v'(r,t) = -rQ(r)E(r,t), (B3) 
2 

where f2(r) is the Keplerian orbital frequency about the binary of mass M given by 



m = J™ (B4) 



Quantity J(r) is the disk angular momentum per unit radius divided by n and is given by 

J(r) = 2r 3 fi(r)S(r). (B5) 

Functions /i(r) and /2( r ) are given by 

/i(r) = 1 P(r)r\ (B6) 
/ 2 (r) = — r 2 + J(r)u7 g (r), (B7) 

where P(r) is the two-dimensional (vertically integrated) disk pressure, S(r) is the disk 
surface density. Quantity w g is the gravitational precession rate of a free particle on an 
eccentric orbit about the equal mass binary, which is given by 

^ (r) = ~^ki d r( r2d ^o(r)) (B8) 

where 

$o(r) = -— X K(x) (B9) 
7rr 
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with 

^W^hr < B10 > 

and K is the complete elliptic integral of the first kind. The terms involving quantities f\ 
and f2 describe the eccentricity propagation and precession, respectively. Quantity 7 is the 
gas adiabatic index. Real function s(r) is the eccentricity injection rate due to some source 
that is assumed to be distributed as a gaussian of width w centered at radius r c and central 
value s c / (y/7rw) 

s(r) = SceXp( - (r " rc)2/w2) . (fill) 
y/irw 

Following along the lines G06 and L10, we adopt the inner boundary condition 

d r E(r h t) = 0. (B12) 

For a Keplerian disk, the divergence of the velocity is proportional to d r E. Consequently, 
this boundary condition is equivalent to requiring that the Lagrangian density perturbation 
near the disk inner edge vanishes. The outer boundary condition is taken to be 

E(r o ,t) = 0. (B13) 

To ensure that the trapping of the eccentricity is not an artifact of the outer boundary 
condition, we test that the resulting mode structure is independent of the outer radius r Q . 

To obtain modes, we write the eccentricity as 

E(r, t) = E(r) exp (iut) (B14) 

where a; is a complex eigenfrequency. The eccentricity distribution equation (Bl) has in- 
finitely many modes and eigenfrequencies. 



C. Eccentricity Distribution Far From Disk Inner Edge 

We determine here the analytic form of the eccentricity distribution far from the inner 
edge of the circumbinary disk, based on the linear theory of Appendix B. In that limit, we 
assume that the local precession rate is small compared with the eigenvalue u defined in 
equation (B14). This assumption is expected to hold, since the gravitational precession rate 
w g declines rapidly with r . We also expect the pressure contribution to the precession rate to 
decline, since the pressure is expected to decline with r. We assume power law dependences 
for the gas sound speed and surface density 

c 2 s (r) = c%x- b \ (CI) 
S(r) = £ a;A (C2) 
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where x = r/a and c s q, S , bi, and b 2 are constants. Equation (Bl) then simplifies to 

x- bl+ V 2 E"(x) + (-& a - b 2 + 3)x~ bl+1/2 E'(x) = QE(x), (C3) 

where we ignored the eccentricity source term that is assumed to be localized near the disk 
inner edge. Dimensionless frequency a; is a constant given by 

7 V C s0 / "bin 

For s>l, equation (C3) has an asymptotic solution in lowest order of the form 

E(x) „ -3^2, m 

where 

* ■ ltI? < C6 > 

c 2 = ^(1 + 260, (C7) 
c 3 = I(9 - 26 1 -46 2 ). (C8) 

o 

In this relation we ignore the overall scale factor for the eccentricity distribution, since we are 
interested in its function form. The appropriate root for y/& is taken such that Re(y/^j) > 0. 

For the case of the simulation, we have that 7 = 1, b\ = (isothermal), c s o = 0.05fibi n o, 
and bj ~ — 0.02fibin*- Quantity 62 = 1 for an isothermal constant a decretion disk, although 
62 is closer to -1 in the outer simulated region during much of the simulation. In any case, 
quantity b 2 has no influence on the dominant factor, the exponential. We adopt b 2 = 1 and 
obtain 

f \ ip/ m exp(-11.3x 1/4 ) 
e(x) = \E(x)\ . (C9) 

We approximate the exponent of equation (C9) by means of a Taylor series about x = 3 
(typical of the simulated outer disk region) and obtain in leading order 

exp (-1.24a;) 

< x ) ^571 ( C1 °) 

and the length scale for the exponential decay is then ~ 0.8a. 
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Fig. 1. — A series of snapshots of the disk surface density at different times. Density contours are 
on a linear scale. The color scale encoding the density (see the color bar for each panel) has twice 
the range in the bottom two panels as in the top two; the white regions are density peaks which 
are off the scale. White dots show the position of the binary; the faint white solid circle shows the 
boundary of the central cut-out; the white dash-dotted circles represent the radii r = 1, 2 and 3a. 
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Fig. 2. — Disk mass profile history. The mass interior to r = 1.0a is given by the solid curve, 1.5a 
the dotted curve, 2.0a the dashed curve, 3a the dash-dotted curve, and 4a the dash-triple-dotted 
curve. 
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Fig. 3. — Time- and shell- integrated surface density (top) and accretion rate (bottom) averaged 
over t = 250-350 (ATi, solid) and t = 350-450 (AT 2 , dashed). The dotted line in the surface 
density plot displays the oc r~ 2 density profile of a 'decretion' disk. 
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Fig. 4. — Top: Time-averaged Reynolds stress (green), Maxwell stress (red) and the total stress 
(black) as functions of radius, solid for ATi and dashed for AT2. Bottom: Same stresses but 
normalized with the local pressure. 
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Fig. 5. — Snapshots of Reynolds stress (left column) and Maxwell stress (right column) at t = 305 
averaged either vertically (top row) or azimuthally (bottom row). The z = ±H and ±2H levels 
are plotted as white dash-dotted lines. 
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Fig. 6. — Time averages (dashed for ATi, solid for AT2) of local (top left), total(top right), and 
specific (bottom left) binary torque as functions of radius. The linear theory prediction using 1 /4 
the time-averaged surface density at the 3 : 2 resonance is shown as a red-dotted curve in the 
upper two panels. The history of the total torque is shown in the bottom right panel. The curve 
is smoothed by boxcar average whose width is twice the binary period. 
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Fig. 7. — Radial derivatives of the different sorts of time-averaged and shell-integrated angular 
momentum flux (dFj/dr) appearing in equation (15) as functions of radius. The averaging period is 
t = 300-320. Net rate of change of the local angular momentum (solid black curve); binary torque 
density (blue dashed curve); Reynolds stress (green dash-dotted curve); Maxwell stress (red dash- 
dotted curve); the sum of Reynolds and Maxwell stresses (cyan dashed curve); angular momentum 
flux due to gas advection (black dashed curve; the inferred flux according to conservation law is 
shown with the magenta dotted curve for comparison). 
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Fig. 8. — Top: Space-time diagram of the local eccentricity (smoothed over T&j n ); color contours 
are logarithmic (see color bar). Bottom: an enlarged section of the space-time diagram (without 
time smoothing) to show short timescale behavior in the inner disk. White regions are eccentricity 
peaks off the scale. 
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Fig. 9. — Time- averaged disk eccentricity as a function of radius at three epochs: t = 150-250 
(black solid curve), t = 250-350 (red dotted curve) and t = 350-450 (blue dashed curve). The gray 
dash-dotted line shows a radial dependence oc exp [2 — 1.3 (r/a)]. 

Table 1. Properties of Simulations of Circumbinary Accretion Disks 



Label 




Type of Simulation 


Resolution 3 


Radial Extent b 


B3D 




MHD 


400 x 160 x 512 


(0.8,16) 


B2D.rin= 


0.8 


Hydrodynamics 


512 x 1024 


(0.8,16) 


B2D.rin= 


1.0 


Hydrodynamics 


472 x 1024 


(1.01,16) 


B2D.rin= 


1.1 


Hydrodynamics 


456 x 1024 


(1.11,16) 


B2D.rin= 


1.2 


Hydrodynamics 


440 x 1024 


(1.22,16) 


B2D.rin= 


1.3 


Hydrodynamics 


428 x 1024 


(1.31,16) 


B2D.rin= 


1.7 


Hydrodynamics 


380 x 1024 


(1.73,16) 



a Resolution is listed as r x 9 x cf) in MHD simulation and r x <p in Hydro- 
dynamic simulations. 

b The azimuthal extent is (0, 2tt) for all simulations. 
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Fig. 10. — Evolution of the disk eccentricity. The dash-dotted lines show linear fits for the 
exponential growth phase and the later saturation phase. 
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Fig. 11. — Phase angle of the disk eccentricity vector as a function of time. 




Fig. 12. — Space-time diagram of the m=l Fourier component of the surface density. Color 
contours are logarithmic (see color bar). 
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Fig. 13. — Three pairs of snapshots of the surface density in the binary's co-rotating frame in 
logarithmic scale (white regions are off the scale) . The density- weighted vertically averaged velocity 
is shown by black arrows. To show how phase-dependent stream effects change secularly over time, 
the time separation between the right and left panels in each pair is a fraction of a binary orbit, 
while the intervals from the first pair to the rest are many orbits. The Roman numerals show the 
quadrants referred in section 3.4. 
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Fig. 14. — Vertically-averaged horizontal magnetic field superposed on surface density contours 
at t = 300.5. Color surface density contours are logarithmic (see color bar), white are peak values 
that is off the scale. 




Fig. 15. — Top: Accretion rate as a function of time (in units of the binary orbital period divided 
by 2ir) during the later stages of the simulation. Middle: Power spectral density of the accretion 
rate as a function of frequency in units of the binary orbital angular frequency £lbm, or 2tt times 
the binary frequency. Bottom: Snapshots of radial mass flux (pv r ) z L z at the times shown by the 
arrows in the top panel. Color contours are in a linear scale (see color bar). Contour lines represent 
the surface density in two groups: black contours show low surface density (10 -4 < £ < 10~ ' 5 ) 
on a logarithmic scale in order to highlight the streams; white contours show high surface density 
(1.4 < £ < 3.0 on a linear scale in order to highlight the lump). 
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Fig. 16. — Disk eccentricity as a function of radius from r = 2a to r = 5a, normalized by its value 
at the inner edge. The heavy solid line is the eccentricity distribution obtained from the simulations 
when the eccentricity evolution is in the linear regime. The middle dashed line is the eccentricity 
distribution determined by linear theory (equation (Bl), with the input parameters appropriate for 
the simulation (see parameter details in § 4.3.1). The thin solid line is the normalized eccentricity 
injection rate distribution, 0.25s(r)/s(r c ). The upper dashed curve is the eccentricity distribution 
for the fundamental free eccentric mode. The lower dashed curve shows the result of a higher 
eccentricity growth rate, about twice as high as in the middle dashed curve. 
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Fig. 17. — The evolution of the disk eccentricity with varying cut-off sizes: r- m = 0.8a (black 
solid curve), 1.0a (blue dash-dotted curve), 1.1a (cyan dash-dot-dotted curve), 1.2a (magenta long 
dashed curve), 1.3a (green dashed curve) and 1.7a (red dotted curve). There is a qualitative change 
when the cut-off becomes larger than ~ 1.2a 




Fig. 19. — Snapshots of the lump concentration and feeding mechanisms at two specific times, 
t = 443 (left) and t = 422 (right), in linear scale (white represents the peak density which is off the 
scale). 



